Create a new analysis directory...
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] "/Users/swvanderlaan/git/CirculatoryHealth/AE_20211201_YAW_SWVANDERLAAN_HDAC9"
[1] "_archived" "1. AE_20211201_YAW_SWVANDERLAAN_HDAC9.nb.html"
[3] "1. AE_20211201_YAW_SWVANDERLAAN_HDAC9.Rmd" "2. SNP_analyses.Rmd"
[5] "20220317.HDAC9.AE_20211201_YAW_SWVANDERLAAN_HDAC9.RData" "20220317.HDAC9.scRNAseq.RData"
[7] "20220318.HDAC9.bulkRNAseq.preparation.RData" "3.1 bulkRNAseq.preparation.nb.html"
[9] "3.1 bulkRNAseq.preparation.Rmd" "3.2 bulkRNAseq.main_analysis.nb.html"
[11] "3.2 bulkRNAseq.main_analysis.Rmd" "3.3 bulkRNAseq.additional_figures.Rmd"
[13] "4. scRNAseq.nb.html" "4. scRNAseq.Rmd"
[15] "AE_20211201_YAW_SWVANDERLAAN_HDAC9.Rproj" "AnalysisPlan"
[17] "HDAC9" "images"
[19] "LICENSE" "README.html"
[21] "README.md" "references.bib"
[23] "renv" "renv.lock"
[25] "scripts" "SNP"
[27] "targets" source(paste0(PROJECT_loc, "/scripts/functions.R"))install.packages.auto("pander")
install.packages.auto("readr")
install.packages.auto("optparse")
install.packages.auto("tools")
install.packages.auto("dplyr")
install.packages.auto("tidyr")
install.packages.auto("naniar")
# To get 'data.table' with 'fwrite' to be able to directly write gzipped-files
# Ref: https://stackoverflow.com/questions/42788401/is-possible-to-use-fwrite-from-data-table-with-gzfile
# install.packages("data.table", repos = "https://Rdatatable.gitlab.io/data.table")
library(data.table)
install.packages.auto("tidyverse")
install.packages.auto("knitr")
install.packages.auto("DT")
install.packages.auto("eeptools")
install.packages.auto("openxlsx")
install.packages.auto("haven")
install.packages.auto("tableone")
install.packages.auto("sjPlot")
install.packages.auto("BlandAltmanLeh")
# Install the devtools package from Hadley Wickham
install.packages.auto('devtools')
# for plotting
install.packages.auto("pheatmap")
install.packages.auto("forestplot")
install.packages.auto("ggplot2")
install.packages.auto("ggpubr")
install.packages.auto("UpSetR")
devtools::install_github("thomasp85/patchwork")Using github PAT from envvar GITHUB_PAT
Skipping install of 'patchwork' from a github remote, the SHA1 (79223d30) has not changed since last install.
Use `force = TRUE` to force installationinstall.packages.auto("GGally")
install.packages.auto("ggcorrplot")
# for Seurat etc
install.packages.auto("GenomicFeatures")
install.packages.auto("GenomicRanges")
install.packages.auto("SummarizedExperiment")
install.packages.auto("DESeq2")
install.packages.auto("org.Hs.eg.db")
install.packages.auto("mygene")
install.packages.auto("TxDb.Hsapiens.UCSC.hg19.knownGene")
install.packages.auto("org.Hs.eg.db")
install.packages.auto("AnnotationDbi")
install.packages.auto("EnsDb.Hsapiens.v86")Loading required package: EnsDb.Hsapiens.v86
Loading required package: ensembldb
Attaching package: 'ensembldb'
The following object is masked from 'package:openxlsx':
addFilter
The following object is masked from 'package:dplyr':
filter
The following object is masked from 'package:stats':
filterinstall.packages.auto("EnhancedVolcano")
Today = format(as.Date(as.POSIXlt(Sys.time())), "%Y%m%d")
Today.Report = format(as.Date(as.POSIXlt(Sys.time())), "%A, %B %d, %Y")
### UtrechtScienceParkColoursScheme
###
### WebsitetoconvertHEXtoRGB:http://hex.colorrrs.com.
### Forsomefunctionsyoushoulddividethesenumbersby255.
###
### No. Color HEX (RGB) CHR MAF/INFO
###---------------------------------------------------------------------------------------
### 1 yellow #FBB820 (251,184,32) => 1 or 1.0>INFO
### 2 gold #F59D10 (245,157,16) => 2
### 3 salmon #E55738 (229,87,56) => 3 or 0.05<MAF<0.2 or 0.4<INFO<0.6
### 4 darkpink #DB003F ((219,0,63) => 4
### 5 lightpink #E35493 (227,84,147) => 5 or 0.8<INFO<1.0
### 6 pink #D5267B (213,38,123) => 6
### 7 hardpink #CC0071 (204,0,113) => 7
### 8 lightpurple #A8448A (168,68,138) => 8
### 9 purple #9A3480 (154,52,128) => 9
### 10 lavendel #8D5B9A (141,91,154) => 10
### 11 bluepurple #705296 (112,82,150) => 11
### 12 purpleblue #686AA9 (104,106,169) => 12
### 13 lightpurpleblue #6173AD (97,115,173/101,120,180) => 13
### 14 seablue #4C81BF (76,129,191) => 14
### 15 skyblue #2F8BC9 (47,139,201) => 15
### 16 azurblue #1290D9 (18,144,217) => 16 or 0.01<MAF<0.05 or 0.2<INFO<0.4
### 17 lightazurblue #1396D8 (19,150,216) => 17
### 18 greenblue #15A6C1 (21,166,193) => 18
### 19 seaweedgreen #5EB17F (94,177,127) => 19
### 20 yellowgreen #86B833 (134,184,51) => 20
### 21 lightmossgreen #C5D220 (197,210,32) => 21
### 22 mossgreen #9FC228 (159,194,40) => 22 or MAF>0.20 or 0.6<INFO<0.8
### 23 lightgreen #78B113 (120,177,19) => 23/X
### 24 green #49A01D (73,160,29) => 24/Y
### 25 grey #595A5C (89,90,92) => 25/XY or MAF<0.01 or 0.0<INFO<0.2
### 26 lightgrey #A2A3A4 (162,163,164) => 26/MT
###
### ADDITIONAL COLORS
### 27 midgrey #D7D8D7
### 28 verylightgrey #ECECEC"
### 29 white #FFFFFF
### 30 black #000000
###----------------------------------------------------------------------------------------------
uithof_color = c("#FBB820","#F59D10","#E55738","#DB003F","#E35493","#D5267B",
"#CC0071","#A8448A","#9A3480","#8D5B9A","#705296","#686AA9",
"#6173AD","#4C81BF","#2F8BC9","#1290D9","#1396D8","#15A6C1",
"#5EB17F","#86B833","#C5D220","#9FC228","#78B113","#49A01D",
"#595A5C","#A2A3A4", "#D7D8D7", "#ECECEC", "#FFFFFF", "#000000")
uithof_color_legend = c("#FBB820", "#F59D10", "#E55738", "#DB003F", "#E35493",
"#D5267B", "#CC0071", "#A8448A", "#9A3480", "#8D5B9A",
"#705296", "#686AA9", "#6173AD", "#4C81BF", "#2F8BC9",
"#1290D9", "#1396D8", "#15A6C1", "#5EB17F", "#86B833",
"#C5D220", "#9FC228", "#78B113", "#49A01D", "#595A5C",
"#A2A3A4", "#D7D8D7", "#ECECEC", "#FFFFFF", "#000000")
### ----------------------------------------------------------------------------Here we perform the main analyses. We first load the prepared RNAseq data, and extract only the relevant genes.
Here we obtain data from the HDAC9 in plaques.
library(openxlsx)
gene_list_df <- read.xlsx(paste0(PROJECT_loc, "/targets/Genes.xlsx"), sheet = "Genes")
target_genes <- unlist(gene_list_df$Gene)
target_genes[1] "HDAC9" "TWIST1" "IL6" "IL1B" We simply load the previously saved RDS-file and extract
the clinical and RNAseq data from that.
AEDB.CEA <- readRDS(file = paste0(OUT_loc, "/20220317.HDAC9.AEDB.CEA.RDS"))
AEDB.CEA$STUDY_NUMBER <- paste0("ae", AEDB.CEA$STUDY_NUMBER)
head(AEDB.CEA$STUDY_NUMBER)[1] "ae1" "ae2" "ae3" "ae4" "ae5" "ae6"AERNASE <- readRDS(file = paste0(OUT_loc, "/20220318.AERNA.CEA.611pts.SE.after_qc.IC_academic.RDS"))AERNASE.clin <- as_tibble(colData(AERNASE))
AERNASE.hdac9 <- subset(AERNASE, subset = (rowRanges(AERNASE)$symbol %in% target_genes))
row.names(AERNASE.hdac9) <- rowData(AERNASE.hdac9)$symbol
temp <- as_tibble(t(assay(AERNASE.hdac9)), rownames = "STUDY_NUMBER")
AERNASE.clin.hdac9 <- as.data.frame(merge(AERNASE.clin, temp, by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = TRUE))
rm(temp)Getting some summary statistics.
cat("Average expression of a random selection of 1000 genes.\n")Average expression of a random selection of 1000 genes.set.seed(141619)
mean(rowMeans(as_tibble(assay(AERNASE)) %>% dplyr::sample_n(1000)))[1] 12.50463cat("\nAverage expression of target genes.\n")
Average expression of target genes.rowMeans(assay(AERNASE.hdac9)) HDAC9 TWIST1 IL1B IL6
28.2062193 5.5859247 5.6072013 0.4157119 cat("\nGene information.\n")
Gene information.rowRanges(AERNASE.hdac9)GRanges object with 4 ranges and 3 metadata columns:
seqnames ranges strand | feature_id symbol GENEBIOTYPE_EnsDb86
<Rle> <IRanges> <Rle> | <character> <character> <character>
HDAC9 7 18086949-19002416 + | ENSG00000048052 HDAC9 protein_coding
TWIST1 7 19020991-19117636 - | ENSG00000122691 TWIST1 protein_coding
IL1B 2 112829751-112836816 - | ENSG00000125538 IL1B protein_coding
IL6 7 22725884-22732002 + | ENSG00000136244 IL6 protein_coding
-------
seqinfo: 338 sequences from an unspecified genome; no seqlengthsAERNASE.hdac9.genedata <- as_tibble(rowRanges(AERNASE.hdac9))
AERNASE.hdac9.genedataThe analyses are focused on three elements:
Age]Gender]Hypertension.composite]DiabetesStatus]SmokerStatus]LDL_final]Med.Statin.LLD]Med.all.antiplatelet]GFR_MDRD]BMI]MedHx_CVD] combination of
[CAD_history, Stroke_history,
Peripheral.interv]stenose]ORdate_year] as we discovered in Van
Lammeren et al. the composition of the plaque and therefore the
Athero-Express Biobank Study has changed over the years. Likely through
changes in lifestyle and primary prevention regimes.We will analyze the data through four different models
In the cross-sectional analysis of plaque MCP1 levels we will focus on the following plaque vulnerability phenotypes:
We inspect the plaque characteristics, and
inverse-rank normal transformation continuous
phenotypes.
# macrophages
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$MAC_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-2.0955 -0.3800 0.3243 0.2763 0.9551 2.6791 15 min_macmean <- min(AERNASE.clin.hdac9$MAC_rankNorm, na.rm = TRUE)
cat(paste0("\nMinimum value % macrophages: ",min_macmean,".\n"))
Minimum value % macrophages: -2.09545532179422.ggpubr::gghistogram(AERNASE.clin.hdac9, "MAC_rankNorm",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "median",
#add_density = TRUE,
rug = TRUE,
#add.params = list(color = "black", linetype = 2),
title = "% macrophages",
xlab = "inverse-rank normalized %",
ggtheme = theme_minimal())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.
Warning: Removed 15 rows containing non-finite values (stat_bin).
# smooth muscle cells
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$SMC_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-2.64518 -0.66839 -0.04095 -0.04085 0.63912 2.55899 17 min_smcmean <- min(AERNASE.clin.hdac9$SMC_rankNorm, na.rm = TRUE)
cat(paste0("\nMinimum value % smooth muscle cells: ",min_smcmean,".\n"))
Minimum value % smooth muscle cells: -2.64517564577702.ggpubr::gghistogram(AERNASE.clin.hdac9, "SMC_rankNorm",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "median",
#add_density = TRUE,
rug = TRUE,
#add.params = list(color = "black", linetype = 2),
title = "% smooth muscle cells",
xlab = "inverse-rank normalized %",
ggtheme = theme_minimal())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.
Warning: Removed 17 rows containing non-finite values (stat_bin).
# vessel density
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$VesselDensity_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-1.9203 -0.4576 0.1109 0.1387 0.7694 3.1703 53 min_vesseldensity <- min(AERNASE.clin.hdac9$VesselDensity_rankNorm, na.rm = TRUE)
min_vesseldensity[1] -1.920282cat(paste0("\nMinimum value number of intraplaque neovessels per 3-4 hotspots: ",min_vesseldensity,".\n"))
Minimum value number of intraplaque neovessels per 3-4 hotspots: -1.92028163343836.ggpubr::gghistogram(AERNASE.clin.hdac9, "VesselDensity_rankNorm",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "median",
#add_density = TRUE,
rug = TRUE,
#add.params = list(color = "black", linetype = 2),
title = "number of intraplaque neovessels per 3-4 hotspots",
xlab = "inverse-rank normalized number",
ggtheme = theme_minimal())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.
Warning: Removed 53 rows containing non-finite values (stat_bin).
Given their strong correlation, we also introduce a macrophages/smooth muscle cell ratio. This is a proxy of the extend to which a plaque is inflammed (‘unstable’) as compared to ‘stable’.
AERNASE.clin.hdac9$MAC_SMC_ratio <- AERNASE.clin.hdac9$MAC_rankNorm / AERNASE.clin.hdac9$SMC_rankNorm
AERNASE.clin.hdac9$MAC_SMC_ratio_rank <- qnorm((rank(AERNASE.clin.hdac9$MAC_SMC_ratio, na.last = "keep") - 0.5) / sum(!is.na(AERNASE.clin.hdac9$MAC_SMC_ratio)))
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$MAC_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-2.0955 -0.3800 0.3243 0.2763 0.9551 2.6791 15 summary(AERNASE.clin.hdac9$SMC_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-2.64518 -0.66839 -0.04095 -0.04085 0.63912 2.55899 17 summary(AERNASE.clin.hdac9$MAC_SMC_ratio_rank) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-3.1410 -0.6732 0.0000 0.0000 0.6732 3.1410 17 ggpubr::gghistogram(AERNASE.clin.hdac9, "MAC_SMC_ratio_rank",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "median",
#add_density = TRUE,
rug = TRUE,
#add.params = list(color = "black", linetype = 2),
title = "macrophages/smooth muscle cells ratio",
xlab = "inverse-rank normalized",
ggtheme = theme_minimal())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.
Warning: Removed 17 rows containing non-finite values (stat_bin).
Warning: call dbDisconnect() when finished working with a connection
library(dplyr)
# calcification
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$Calc.bin) no/minor moderate/heavy NA's
283 320 8 contrasts(AERNASE.clin.hdac9$Calc.bin) moderate/heavy
no/minor 0
moderate/heavy 1AERNASE.clin.hdac9$CalcificationPlaque <- as.factor(AERNASE.clin.hdac9$Calc.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(CalcificationPlaque)) %>%
dplyr::group_by(Gender, CalcificationPlaque) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "CalcificationPlaque", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Calcification",
xlab = "calcification",
ggtheme = theme_minimal())
rm(df)
# collagen
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$Collagen.bin) no/minor moderate/heavy NA's
117 485 9 contrasts(AERNASE.clin.hdac9$Collagen.bin) moderate/heavy
no/minor 0
moderate/heavy 1AERNASE.clin.hdac9$CollagenPlaque <- as.factor(AERNASE.clin.hdac9$Collagen.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(CollagenPlaque)) %>%
dplyr::group_by(Gender, CollagenPlaque) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "CollagenPlaque", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Collagen",
xlab = "collagen",
ggtheme = theme_minimal())
rm(df)
# fat 10%
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$Fat.bin_10) <10% >10% NA's
143 460 8 contrasts(AERNASE.clin.hdac9$Fat.bin_10) >10%
<10% 0
>10% 1AERNASE.clin.hdac9$Fat10Perc <- as.factor(AERNASE.clin.hdac9$Fat.bin_10)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(Fat10Perc)) %>%
dplyr::group_by(Gender, Fat10Perc) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "Fat10Perc", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Intraplaque fat",
xlab = "intraplaque fat",
ggtheme = theme_minimal())
rm(df)
# macrophages binned
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$Macrophages.bin) no/minor moderate/heavy NA's
261 339 11 contrasts(AERNASE.clin.hdac9$Macrophages.bin) moderate/heavy
no/minor 0
moderate/heavy 1AERNASE.clin.hdac9$MAC_binned <- as.factor(AERNASE.clin.hdac9$Macrophages.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(MAC_binned)) %>%
dplyr::group_by(Gender, MAC_binned) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "MAC_binned", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Macrophages (binned)",
xlab = "Macrophages",
ggtheme = theme_minimal())
rm(df)
# SMC binned
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$SMC.bin) no/minor moderate/heavy NA's
193 407 11 contrasts(AERNASE.clin.hdac9$SMC.bin) moderate/heavy
no/minor 0
moderate/heavy 1AERNASE.clin.hdac9$SMC_binned <- as.factor(AERNASE.clin.hdac9$SMC.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(SMC_binned)) %>%
dplyr::group_by(Gender, SMC_binned) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "SMC_binned", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "SMC (binned)",
xlab = "SMC",
ggtheme = theme_minimal())
rm(df)
# IPH
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$IPH.bin) no yes NA's
230 372 9 contrasts(AERNASE.clin.hdac9$IPH.bin) yes
no 0
yes 1AERNASE.clin.hdac9$IPH <- as.factor(AERNASE.clin.hdac9$IPH.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(IPH)) %>%
dplyr::group_by(Gender, IPH) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "IPH", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Intraplaque hemorrhage",
xlab = "intraplaque hemorrhage",
ggtheme = theme_minimal())
rm(df)
# Symptoms
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$AsymptSympt)Ocular and others Symptomatic
184 427 contrasts(AERNASE.clin.hdac9$AsymptSympt) Symptomatic
Ocular and others 0
Symptomatic 1AERNASE.clin.hdac9$AsymptSympt <- as.factor(AERNASE.clin.hdac9$AsymptSympt)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(AsymptSympt)) %>%
dplyr::group_by(Gender, AsymptSympt) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "AsymptSympt", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Symptoms",
xlab = "symptoms",
ggtheme = theme_minimal())
rm(df)Here we compare the HDAC9 plaque gene expression levels.
p1 <- ggpubr::ggscatter(AERNASE.clin.hdac9,
x = "ORyear",
y = "HDAC9",
color = "#1290D9",
# fill = "Gender",
# palette = c("#1290D9", "#DB003F"),
add = "reg.line",
add.params = list(color = "black", linetype = 2),
cor.coef = TRUE, cor.method = "spearman",
xlab = "year of surgery",
ylab = "experiment 1",
title = "HDAC9 (normalized expression)",
ggtheme = theme_minimal())
p1 `geom_smooth()` using formula 'y ~ x'
rm(p1)In this section we make some variables to assist with analysis.
We fix the diabetes status variable.
# Fix diabetes
attach(AEDB.CEA)
AEDB.CEA[,"DiabetesStatus"] <- NA
AEDB.CEA$DiabetesStatus[DM.composite == -999] <- NA
AEDB.CEA$DiabetesStatus[DM.composite == 1] <- "Control (no Diabetes Dx/Med)"
AEDB.CEA$DiabetesStatus[DM.composite == 2] <- "Diabetes"
detach(AEDB.CEA)
table(AEDB.CEA$DM.composite, AEDB.CEA$DiabetesStatus)
Control (no Diabetes Dx/Med) Diabetes
1 1900 0
2 0 596# AEDB.CEA.temp <- subset(AEDB.CEA, select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "DM.composite", "DiabetesStatus"))
# require(labelled)
# AEDB.CEA.temp$Gender <- to_factor(AEDB.CEA.temp$Gender)
# AEDB.CEA.temp$Hospital <- to_factor(AEDB.CEA.temp$Hospital)
# AEDB.CEA.temp$Artery_summary <- to_factor(AEDB.CEA.temp$Artery_summary)
# AEDB.CEA.temp$DiabetesStatus <- to_factor(AEDB.CEA.temp$DiabetesStatus)
#
# DT::datatable(AEDB.CEA.temp[1:10,], caption = "Excerpt of the whole AEDB.CEA.", rownames = FALSE)
#
# rm(AEDB.CEA.temp)We will also fix a history of CAD, stroke or peripheral intervention
status variable. This will be based on CAD_history,
Stroke_history, and Peripheral.interv
# AEDB.CEA$CAD_history
# AEDB.CEA$Stroke_history
# AEDB.CEA$Peripheral.interv
# Fix diabetes
attach(AEDB.CEA)
AEDB.CEA[,"MedHx_CVD"] <- NA
AEDB.CEA$MedHx_CVD[CAD_history == "No history CAD" | Stroke_history == 1 | Peripheral.interv == "no"] <- "No"
AEDB.CEA$MedHx_CVD[CAD_history == "History CAD" | Stroke_history == 2 | Peripheral.interv == "yes"] <- "yes"
detach(AEDB.CEA)
table(AEDB.CEA$CAD_history)
No history CAD History CAD
1703 791 table(AEDB.CEA$Stroke_history)
1 2
1691 821 table(AEDB.CEA$Peripheral.interv)
no yes
1971 515 table(AEDB.CEA$MedHx_CVD)
No yes
909 1603 # AEDB.CEA.temp <- subset(AEDB.CEA, select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "diet810", "AlcoholUse"))
# require(labelled)
# AEDB.CEA.temp$Gender <- to_factor(AEDB.CEA.temp$Gender)
# AEDB.CEA.temp$Hospital <- to_factor(AEDB.CEA.temp$Hospital)
# AEDB.CEA.temp$Artery_summary <- to_factor(AEDB.CEA.temp$Artery_summary)
# AEDB.CEA.temp$AlcoholUse <- to_factor(AEDB.CEA.temp$AlcoholUse)
#
# DT::datatable(AEDB.CEA.temp[1:10,], caption = "Excerpt of the whole AEDB.CEA.", rownames = FALSE)
#
# rm(AEDB.CEA.temp)AERNASE.clin.hdac9$DiabetesStatus <- NULL
AERNASE.clin.hdac9$MedHx_CVD <- NULL
AERNASE.clin.hdac9 <- merge(AERNASE.clin.hdac9,
subset(AEDB.CEA, select = c("STUDY_NUMBER", "dateok", "DiabetesStatus", "MedHx_CVD")),
by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = TRUE)
AERNASE.clin.hdac9$MedHx_CVD <- as_factor(AERNASE.clin.hdac9$MedHx_CVD)
AERNASE.clin.hdac9$DiabetesStatus <- as_factor(AERNASE.clin.hdac9$DiabetesStatus)
study.samplesize = nrow(AERNASE.clin.hdac9) # study.samplesize is an expected variable in the GLM.BIN() GLM.CON() functions
TRAITS.TARGET.RANK = c("HDAC9")
TRAITS.CON.RANK = c("MAC_rankNorm", "SMC_rankNorm", "MAC_SMC_ratio_rank", "VesselDensity_rankNorm")
TRAITS.BIN = c("CalcificationPlaque", "CollagenPlaque", "Fat10Perc", "IPH",
"MAC_binned", "SMC_binned")
# "Hospital",
# "Age", "Gender",
# "TC_final", "LDL_final", "HDL_final", "TG_final",
# "systolic", "diastoli", "GFR_MDRD", "BMI",
# "KDOQI", "BMI_WHO",
# "SmokerCurrent", "eCigarettes", "ePackYearsSmoking",
# "DiabetesStatus", "Hypertension.composite",
# "Hypertension.drugs", "Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
# "Stroke_Dx", "sympt", "Symptoms.5G", "restenos",
# "EP_composite", "EP_composite_time",
# "macmean0", "smcmean0", "Macrophages.bin", "SMC.bin",
# "neutrophils", "Mast_cells_plaque",
# "IPH.bin", "vessel_density_averaged",
# "Calc.bin", "Collagen.bin",
# "Fat.bin_10", "Fat.bin_40", "OverallPlaquePhenotype",
# "IL6_pg_ug_2015", "MCP1_pg_ug_2015",
# "QC2018_FILTER", "CHIP", "SAMPLE_TYPE",
# "CAD_history", "Stroke_history", "Peripheral.interv",
# "stenose"
# 1. Age (continuous in 1-year increment). [Age]
# 2. Sex (male vs. female). [Gender]
# 3. Presence of hypertension at baseline (defined either as history of hypertension, SBP ≥140 mm Hg, DBP ≥90 mm Hg, or prescription of antihypertensive medications). [Hypertension.composite]
# 4. Presence of diabetes mellitus at baseline (defined either as a history of diabetes, administration of glucose lowering medication, HbA1c ≥6.5%, fasting glucose ≥126 mg/dl, .or random glucose levels ≥200 mg/dl). [DiabetesStatus]
# 5. Smoking (current, ex-, never). [SmokerCurrent]
# 6. LDL-C levels (continuous). [LDL_final]
# 7. Use of lipid-lowering drugs. [Med.Statin.LLD]
# 8. Use of antiplatelet drugs. [Med.all.antiplatelet]
# 9. eGFR (continuous). [GFR_MDRD]
# 10. BMI (continuous). [BMI]
# 11. History of cardiovascular disease (stroke, coronary artery disease, peripheral artery disease). [MedHx_CVD] combinatino of: [CAD_history, Stroke_history, Peripheral.interv]
# 12. Level of stenosis (50-70% vs. 70-99%). [stenose]
# Models
# Model 1: adjusted for age and sex
# Model 2: adjusted for age, sex, hypertension, diabetes, smoking, LDL-C levels, lipid-lowering drugs, antiplatelet drugs, eGFR, BMI, history of CVD, level of stenosis,
AERNASE.clin.hdac9$ORdate_epoch <- as.numeric(AERNASE.clin.hdac9$dateok)
AERNASE.clin.hdac9$ORdate_year <- AERNASE.clin.hdac9$ORyear
cat("Summary of 'year of surgery' in 'epoch' (); coded as `numeric()`\n")Summary of 'year of surgery' in 'epoch' (); coded as `numeric()`summary(AERNASE.clin.hdac9$ORdate_epoch) Min. 1st Qu. Median Mean 3rd Qu. Max.
1.324e+10 1.332e+10 1.338e+10 1.339e+10 1.346e+10 1.367e+10 cat("\nSummary of 'year of surgery' in 'years' (); coded as `factor()`\n")
Summary of 'year of surgery' in 'years' (); coded as `factor()`table(AERNASE.clin.hdac9$ORdate_year)
2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022
32 61 65 80 81 64 62 67 34 31 22 5 3 3 1 0 0 0 0 0 0 # COVARIATES_M1 = c("Age", "Gender", "ORdate_year")
COVARIATES_M1 = c("Age", "Gender", "ORdate_epoch")
COVARIATES_M2 = c(COVARIATES_M1,
"Hypertension.composite", "DiabetesStatus",
"SmokerStatus",
# "SmokerCurrent",
"Med.Statin.LLD", "Med.all.antiplatelet",
"GFR_MDRD", "BMI",
# "CAD_history", "Stroke_history", "Peripheral.interv",
"MedHx_CVD",
"stenose")
str(AERNASE.clin.hdac9)'data.frame': 611 obs. of 74 variables:
$ STUDY_NUMBER : chr "ae1" "ae1026" "ae1029" "ae1032" ...
$ SampleType : chr "plaque" "plaque" "plaque" "plaque" ...
$ RNAseqType : chr "3' RNAseq" "3' RNAseq" "3' RNAseq" "3' RNAseq" ...
$ Hospital : Factor w/ 2 levels "St. Antonius, Nieuwegein",..: 1 2 2 2 2 2 2 1 2 2 ...
$ ORyear : Factor w/ 21 levels "2002","2003",..: 1 4 4 4 4 4 5 2 4 6 ...
$ Artery_summary : Factor w/ 9 levels "No artery known (yet), no surgery (patient ill, died, exited study), re-numbered to AAA",..: 2 2 2 2 2 2 2 2 2 2 ...
$ Age : num 29 46 23 15 41 37 43 44 26 25 ...
$ Gender : Factor w/ 2 levels "female","male": 2 2 2 1 1 1 2 2 2 1 ...
$ TC_final : num 4.48 NA 4.2 NA 2.92 ...
$ LDL_final : num 2.94 NA NA NA 0.65 ...
$ HDL_final : num 1 NA NA NA 0.87 ...
$ TG_final : num 1.21 NA NA NA 3.07 ...
$ systolic : num 129 71 61 81 15 91 71 76 51 61 ...
$ diastoli : num 51 41 44 51 16 51 51 31 21 13 ...
$ GFR_MDRD : num 491 951 1499 2580 135 ...
$ BMI : num 879 682 208 245 1048 ...
$ KDOQI : Factor w/ 5 levels "Normal kidney function",..: 3 2 2 1 3 4 2 3 3 3 ...
$ BMI_WHO : Factor w/ 4 levels "Underweight",..: 3 3 2 2 4 3 3 2 3 2 ...
$ SmokerStatus : Factor w/ 3 levels "Current smoker",..: 2 NA 1 1 3 1 3 3 3 1 ...
$ AlcoholUse : Factor w/ 1 level "Yes": 1 NA 1 NA NA 1 NA 1 NA 1 ...
$ Hypertension.selfreport : Factor w/ 2 levels "no","yes": 2 1 1 1 2 2 2 2 2 2 ...
$ Hypertension.selfreportdrug: Factor w/ 2 levels "no","yes": 2 1 1 1 2 1 2 2 2 2 ...
$ Hypertension.composite : Factor w/ 2 levels "no","yes": 2 2 1 1 2 2 2 2 2 2 ...
$ Hypertension.drugs : Factor w/ 2 levels "no","yes": 1 2 1 1 2 2 2 2 2 2 ...
$ Med.anticoagulants : Factor w/ 2 levels "no","yes": 1 2 1 1 1 2 1 2 1 1 ...
$ Med.all.antiplatelet : Factor w/ 2 levels "no","yes": 2 1 2 2 2 1 2 2 2 2 ...
$ Med.Statin.LLD : Factor w/ 2 levels "no","yes": 1 2 2 1 2 2 2 1 2 2 ...
$ Stroke_Dx : Factor w/ 2 levels "No stroke diagnosed",..: 1 1 1 1 1 1 1 1 1 1 ...
$ sympt : Factor w/ 17 levels "Asymptomatic",..: 2 2 4 2 5 3 1 1 1 1 ...
$ Symptoms.5G : Factor w/ 5 levels "Ocular","Other",..: 4 4 1 4 2 4 5 5 5 5 ...
$ AsymptSympt : Factor w/ 2 levels "Ocular and others",..: 2 2 1 2 1 2 2 2 2 2 ...
$ AsymptSympt2G : Factor w/ 1 level "Symptomatic": 1 1 1 1 1 1 1 1 1 1 ...
$ Symptoms.Update2G : Factor w/ 1 level "Symptomatic": 1 1 1 1 1 1 1 1 1 1 ...
$ Symptoms.Update3G : Factor w/ 1 level "Symptomatic": 1 1 1 1 1 1 1 1 1 1 ...
$ restenos : Factor w/ 3 levels "de novo","restenosis",..: 1 1 1 1 1 1 1 1 1 1 ...
$ stenose : Factor w/ 9 levels "0-49%","50-70%",..: 4 2 3 3 4 3 3 4 3 3 ...
$ CAD_history : Factor w/ 2 levels "No history CAD",..: 2 1 1 1 1 1 1 1 1 1 ...
$ PAOD : Factor w/ 2 levels "no","yes": 2 2 1 1 1 1 1 1 1 2 ...
$ Peripheral.interv : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 2 ...
$ EP_composite : Factor w/ 2 levels "No composite endpoints",..: 1 2 2 1 2 2 2 1 1 1 ...
$ EP_composite_time : num 12.258 0.208 0.685 1.288 2.603 ...
$ EP_major : Factor w/ 2 levels "No major events (endpoints)",..: 1 1 1 1 2 2 2 1 1 1 ...
$ EP_major_time : num 12.258 0.647 2.989 1.288 2.603 ...
$ MAC_rankNorm : num -0.0697 -0.7511 0.9114 -0.38 -1.0968 ...
$ SMC_rankNorm : num 1.522 0.58 0.673 1.256 0.342 ...
$ Macrophages.bin : Factor w/ 2 levels "no/minor","moderate/heavy": 2 1 2 1 1 2 2 1 2 1 ...
$ SMC.bin : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 2 2 2 1 2 2 2 ...
$ Neutrophils_rankNorm : num NA 0.591 0.118 -1.618 -1.618 ...
$ MastCells_rankNorm : num NA 0.633 NA NA NA ...
$ IPH.bin : Factor w/ 2 levels "no","yes": 1 2 2 1 1 1 1 2 1 1 ...
$ VesselDensity_rankNorm : num -0.723 0.769 0.854 1.14 -0.861 ...
$ Calc.bin : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 1 1 1 2 2 2 1 ...
$ Collagen.bin : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 2 2 2 1 2 2 2 ...
$ Fat.bin_10 : Factor w/ 2 levels " <10%"," >10%": 2 2 2 2 1 1 2 1 2 1 ...
$ Fat.bin_40 : Factor w/ 2 levels "<40%",">40%": 1 1 1 1 1 1 2 1 1 1 ...
$ OverallPlaquePhenotype : Factor w/ 3 levels "atheromatous",..: 1 2 1 1 2 2 NA 2 1 2 ...
$ Plaque_Vulnerability_Index : Factor w/ 6 levels "0","1","2","3",..: 3 3 4 2 1 2 5 2 3 1 ...
$ HDAC9 : int 24 11 7 5 13 46 14 0 30 12 ...
$ TWIST1 : int 1 13 5 0 1 1 3 1 19 2 ...
$ IL1B : int 2 4 2 1 1 2 5 2 43 2 ...
$ IL6 : int 0 1 0 0 0 0 0 0 5 0 ...
$ MAC_SMC_ratio : num -0.0458 -1.2947 1.3537 -0.3026 -3.2041 ...
$ MAC_SMC_ratio_rank : num -0.215 -0.88 0.728 -0.36 -1.31 ...
$ CalcificationPlaque : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 1 1 1 2 2 2 1 ...
$ CollagenPlaque : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 2 2 2 1 2 2 2 ...
$ Fat10Perc : Factor w/ 2 levels " <10%"," >10%": 2 2 2 2 1 1 2 1 2 1 ...
$ MAC_binned : Factor w/ 2 levels "no/minor","moderate/heavy": 2 1 2 1 1 2 2 1 2 1 ...
$ SMC_binned : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 2 2 2 1 2 2 2 ...
$ IPH : Factor w/ 2 levels "no","yes": 1 2 2 1 1 1 1 2 1 1 ...
$ dateok : num 1.32e+10 1.33e+10 1.33e+10 1.33e+10 1.33e+10 ...
$ DiabetesStatus : Factor w/ 2 levels "Control (no Diabetes Dx/Med)",..: 1 1 1 1 2 1 1 1 1 1 ...
$ MedHx_CVD : Factor w/ 2 levels "yes","No": 1 2 1 2 2 1 2 2 2 1 ...
$ ORdate_epoch : num 1.32e+10 1.33e+10 1.33e+10 1.33e+10 1.33e+10 ...
$ ORdate_year : Factor w/ 21 levels "2002","2003",..: 1 4 4 4 4 4 5 2 4 6 ...# COVARIATES_M3 = c(COVARIATES_M2, "LDL_final")
# COVARIATES_M4 = c(COVARIATES_M2, "hsCRP_plasma")In this model we correct for Age, Gender, and year of surgery.
Here we use the inverse-rank normalized data - visually this is more normally distributed.
detach("package:EnsDb.Hsapiens.v86", unload = TRUE)
detach("package:ensembldb", unload = TRUE)Analysis of continuous/quantitative plaque traits as a function of plaque HDAC9 expression levels.
library(dplyr)
library(MASS)
GLM.results <- data.frame(matrix(NA, ncol = 15, nrow = 0))
cat("Running linear regression...\n")Running linear regression...for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
for (trait in 1:length(TRAITS.CON.RANK)) {
TRAIT = TRAITS.CON.RANK[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M1) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# # for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
### univariate
# fit <- lm(currentDF[,TARGET] ~ currentDF[,TRAIT] + Age + Gender + ORdate_year, data = currentDF)
fit <- lm(currentDF[,TARGET] ~ currentDF[,TRAIT] + Age + Gender + ORdate_epoch, data = currentDF)
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 15, nrow = 0))
GLM.results.TEMP[1,] = GLM.CON(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
Analysis of HDAC9.
- processing MAC_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ 1, data = currentDF)
Coefficients:
(Intercept)
28.22
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-29.01 -17.56 -7.40 8.71 420.59
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.289e+02 2.147e+02 -0.600 0.549
currentDF[, TRAIT] 2.065e-01 1.399e+00 0.148 0.883
Age 1.415e-01 1.500e-01 0.943 0.346
Gendermale -1.241e+00 3.096e+00 -0.401 0.689
ORdate_epoch 1.144e-08 1.607e-08 0.712 0.477
Residual standard error: 32.39 on 591 degrees of freedom
Multiple R-squared: 0.002791, Adjusted R-squared: -0.003959
F-statistic: 0.4135 on 4 and 591 DF, p-value: 0.799
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_rankNorm
Effect size...............: 0.206547
Standard error............: 1.398569
Odds ratio (effect size)..: 1.229
Lower 95% CI..............: 0.079
Upper 95% CI..............: 19.063
T-value...................: 0.147684
P-value...................: 0.8826422
R^2.......................: 0.002791
Adjusted r^2..............: -0.003959
Sample size of AE DB......: 611
Sample size of model......: 596
Missing data %............: 2.454992
- processing SMC_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT], data = currentDF)
Coefficients:
(Intercept) currentDF[, TRAIT]
27.975 -4.619
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-34.98 -17.20 -7.50 7.82 420.01
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.827e+01 2.148e+02 -0.132 0.89536
currentDF[, TRAIT] -4.555e+00 1.471e+00 -3.097 0.00205 **
Age 6.111e-02 1.508e-01 0.405 0.68556
Gendermale -2.257e+00 3.090e+00 -0.730 0.46549
ORdate_epoch 4.172e-09 1.608e-08 0.260 0.79533
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 32.16 on 589 degrees of freedom
Multiple R-squared: 0.01891, Adjusted R-squared: 0.01225
F-statistic: 2.838 on 4 and 589 DF, p-value: 0.02377
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' SMC_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: SMC_rankNorm
Effect size...............: -4.554888
Standard error............: 1.470541
Odds ratio (effect size)..: 0.011
Lower 95% CI..............: 0.001
Upper 95% CI..............: 0.188
T-value...................: -3.097424
P-value...................: 0.002045031
R^2.......................: 0.018908
Adjusted r^2..............: 0.012245
Sample size of AE DB......: 611
Sample size of model......: 594
Missing data %............: 2.782324
- processing MAC_SMC_ratio_rank
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT], data = currentDF)
Coefficients:
(Intercept) currentDF[, TRAIT]
28.163 -2.839
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-33.32 -17.33 -7.61 7.49 418.30
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.661e+02 2.131e+02 -0.779 0.4361
currentDF[, TRAIT] -2.929e+00 1.330e+00 -2.202 0.0281 *
Age 1.336e-01 1.494e-01 0.895 0.3713
Gendermale -1.655e+00 3.091e+00 -0.535 0.5926
ORdate_epoch 1.426e-08 1.596e-08 0.893 0.3720
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 32.29 on 589 degrees of freedom
Multiple R-squared: 0.01107, Adjusted R-squared: 0.004349
F-statistic: 1.648 on 4 and 589 DF, p-value: 0.1607
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_SMC_ratio_rank ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_SMC_ratio_rank
Effect size...............: -2.928947
Standard error............: 1.33041
Odds ratio (effect size)..: 0.053
Lower 95% CI..............: 0.004
Upper 95% CI..............: 0.725
T-value...................: -2.201537
P-value...................: 0.02808507
R^2.......................: 0.011065
Adjusted r^2..............: 0.004349
Sample size of AE DB......: 611
Sample size of model......: 594
Missing data %............: 2.782324
- processing VesselDensity_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ 1, data = currentDF)
Coefficients:
(Intercept)
28.58
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-29.63 -17.87 -7.71 8.54 420.16
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.649e+02 2.272e+02 -0.726 0.468
currentDF[, TRAIT] -5.828e-01 1.485e+00 -0.392 0.695
Age 1.461e-01 1.583e-01 0.923 0.356
Gendermale -1.405e+00 3.260e+00 -0.431 0.667
ORdate_epoch 1.417e-08 1.701e-08 0.833 0.405
Residual standard error: 33.08 on 553 degrees of freedom
Multiple R-squared: 0.003666, Adjusted R-squared: -0.003541
F-statistic: 0.5087 on 4 and 553 DF, p-value: 0.7294
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' VesselDensity_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: VesselDensity_rankNorm
Effect size...............: -0.582804
Standard error............: 1.485197
Odds ratio (effect size)..: 0.558
Lower 95% CI..............: 0.03
Upper 95% CI..............: 10.259
T-value...................: -0.392408
P-value...................: 0.6949077
R^2.......................: 0.003666
Adjusted r^2..............: -0.003541
Sample size of AE DB......: 611
Sample size of model......: 558
Missing data %............: 8.674304 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"T-value", "P-value", "r^2", "r^2_adj", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`T-value` <- as.numeric(GLM.results$`T-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2` <- as.numeric(GLM.results$`r^2`)
GLM.results$`r^2_adj` <- as.numeric(GLM.results$`r^2_adj`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
library(openxlsx)
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Con.Uni.Protein.PlaquePhenotypes.RANK.MODEL1.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Con.Uni.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)Analysis of binary plaque traits as a function of plaque HDAC9 expression levels.
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
for (trait in 1:length(TRAITS.BIN)) {
TRAIT = TRAITS.BIN[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M1) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_year,
# data = currentDF, family = binomial(link = "logit"))
fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_epoch,
data = currentDF, family = binomial(link = "logit"))
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 16, nrow = 0))
GLM.results.TEMP[1,] = GLM.BIN(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
Analysis of HDAC9.
- processing CalcificationPlaque
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale ORdate_epoch
9.089e+01 5.593e-03 -3.617e-01 -6.770e-09
Degrees of Freedom: 602 Total (i.e. Null); 599 Residual
Null Deviance: 833.7
Residual Deviance: 780.8 AIC: 788.8
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.856 -1.115 0.768 1.052 1.710
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 9.127e+01 1.403e+01 6.506 7.71e-11 ***
currentDF[, TARGET] 5.549e-03 3.146e-03 1.764 0.0778 .
Age 4.413e-03 9.607e-03 0.459 0.6460
Gendermale -3.612e-01 2.012e-01 -1.796 0.0725 .
ORdate_epoch -6.811e-09 1.050e-09 -6.485 8.88e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 833.66 on 602 degrees of freedom
Residual deviance: 780.57 on 598 degrees of freedom
AIC: 790.57
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' CalcificationPlaque ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: CalcificationPlaque
Effect size...............: 0.005549
Standard error............: 0.003146
Odds ratio (effect size)..: 1.006
Lower 95% CI..............: 0.999
Upper 95% CI..............: 1.012
Z-value...................: 1.763675
P-value...................: 0.07778677
Hosmer and Lemeshow r^2...: 0.063691
Cox and Snell r^2.........: 0.084289
Nagelkerke's pseudo r^2...: 0.112527
Sample size of AE DB......: 611
Sample size of model......: 603
Missing data %............: 1.309329
- processing CollagenPlaque
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ 1, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept)
1.422
Degrees of Freedom: 601 Total (i.e. Null); 601 Residual
Null Deviance: 592.9
Residual Deviance: 592.9 AIC: 594.9
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9490 0.5987 0.6368 0.6722 1.1578
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.634e+01 1.615e+01 1.012 0.312
currentDF[, TARGET] -3.491e-03 2.812e-03 -1.242 0.214
Age 4.386e-03 1.162e-02 0.378 0.706
Gendermale 1.249e-01 2.380e-01 0.525 0.600
ORdate_epoch -1.125e-09 1.208e-09 -0.931 0.352
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 592.94 on 601 degrees of freedom
Residual deviance: 590.28 on 597 degrees of freedom
AIC: 600.28
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' CollagenPlaque ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: CollagenPlaque
Effect size...............: -0.003491
Standard error............: 0.002812
Odds ratio (effect size)..: 0.997
Lower 95% CI..............: 0.991
Upper 95% CI..............: 1.002
Z-value...................: -1.241602
P-value...................: 0.2143835
Hosmer and Lemeshow r^2...: 0.004473
Cox and Snell r^2.........: 0.004396
Nagelkerke's pseudo r^2...: 0.007016
Sample size of AE DB......: 611
Sample size of model......: 602
Missing data %............: 1.472995
- processing Fat10Perc
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale ORdate_epoch
8.888e+01 1.881e-02 7.717e-01 -6.623e-09
Degrees of Freedom: 602 Total (i.e. Null); 599 Residual
Null Deviance: 660.6
Residual Deviance: 607.9 AIC: 615.9
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.5116 0.2827 0.5812 0.7565 1.3154
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.998e+01 1.609e+01 5.593 2.24e-08 ***
currentDF[, TARGET] 1.866e-02 5.470e-03 3.412 0.000644 ***
Age 1.226e-02 1.109e-02 1.106 0.268904
Gendermale 7.758e-01 2.237e-01 3.468 0.000524 ***
ORdate_epoch -6.736e-09 1.205e-09 -5.592 2.25e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 660.61 on 602 degrees of freedom
Residual deviance: 606.67 on 598 degrees of freedom
AIC: 616.67
Number of Fisher Scoring iterations: 5
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' Fat10Perc ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: Fat10Perc
Effect size...............: 0.018664
Standard error............: 0.00547
Odds ratio (effect size)..: 1.019
Lower 95% CI..............: 1.008
Upper 95% CI..............: 1.03
Z-value...................: 3.412234
P-value...................: 0.0006443268
Hosmer and Lemeshow r^2...: 0.081656
Cox and Snell r^2.........: 0.085572
Nagelkerke's pseudo r^2...: 0.128556
Sample size of AE DB......: 611
Sample size of model......: 603
Missing data %............: 1.309329
- processing IPH
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender + ORdate_epoch,
family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) Gendermale ORdate_epoch
8.801e+01 6.163e-01 -6.570e-09
Degrees of Freedom: 601 Total (i.e. Null); 599 Residual
Null Deviance: 800.7
Residual Deviance: 752.1 AIC: 758.1
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9201 -1.1952 0.7258 0.9514 1.7063
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.891e+01 1.419e+01 6.267 3.69e-10 ***
currentDF[, TARGET] 7.551e-04 2.829e-03 0.267 0.78956
Age 8.173e-03 9.843e-03 0.830 0.40631
Gendermale 6.192e-01 2.010e-01 3.081 0.00207 **
ORdate_epoch -6.660e-09 1.063e-09 -6.267 3.67e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 800.74 on 601 degrees of freedom
Residual deviance: 751.31 on 597 degrees of freedom
AIC: 761.31
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' IPH ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: IPH
Effect size...............: 0.000755
Standard error............: 0.002829
Odds ratio (effect size)..: 1.001
Lower 95% CI..............: 0.995
Upper 95% CI..............: 1.006
Z-value...................: 0.266882
P-value...................: 0.7895597
Hosmer and Lemeshow r^2...: 0.061732
Cox and Snell r^2.........: 0.078831
Nagelkerke's pseudo r^2...: 0.107171
Sample size of AE DB......: 611
Sample size of model......: 602
Missing data %............: 1.472995
- processing MAC_binned
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) Gendermale
-0.09663 0.47469
Degrees of Freedom: 599 Total (i.e. Null); 598 Residual
Null Deviance: 821.6
Residual Deviance: 815.5 AIC: 819.5
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7639 -1.2859 0.9725 1.0439 1.3126
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.197e+01 1.303e+01 0.919 0.3582
currentDF[, TARGET] 3.936e-03 2.974e-03 1.324 0.1856
Age 1.141e-03 9.342e-03 0.122 0.9028
Gendermale 4.921e-01 1.926e-01 2.555 0.0106 *
ORdate_epoch -9.135e-10 9.757e-10 -0.936 0.3491
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 821.61 on 599 degrees of freedom
Residual deviance: 812.70 on 595 degrees of freedom
AIC: 822.7
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_binned ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_binned
Effect size...............: 0.003936
Standard error............: 0.002974
Odds ratio (effect size)..: 1.004
Lower 95% CI..............: 0.998
Upper 95% CI..............: 1.01
Z-value...................: 1.323626
P-value...................: 0.1856274
Hosmer and Lemeshow r^2...: 0.010844
Cox and Snell r^2.........: 0.01474
Nagelkerke's pseudo r^2...: 0.019766
Sample size of AE DB......: 611
Sample size of model......: 600
Missing data %............: 1.800327
- processing SMC_binned
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age, family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Age
1.677785 -0.004282 -0.023224
Degrees of Freedom: 599 Total (i.e. Null); 597 Residual
Null Deviance: 753.7
Residual Deviance: 745.4 AIC: 751.4
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8005 -1.3987 0.8089 0.8977 1.7277
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.695e+00 1.383e+01 0.412 0.6805
currentDF[, TARGET] -4.322e-03 2.741e-03 -1.577 0.1149
Age -2.301e-02 1.013e-02 -2.272 0.0231 *
Gendermale -2.441e-01 2.111e-01 -1.156 0.2475
ORdate_epoch -2.865e-10 1.036e-09 -0.277 0.7821
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 753.74 on 599 degrees of freedom
Residual deviance: 743.90 on 595 degrees of freedom
AIC: 753.9
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' SMC_binned ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: SMC_binned
Effect size...............: -0.004322
Standard error............: 0.002741
Odds ratio (effect size)..: 0.996
Lower 95% CI..............: 0.99
Upper 95% CI..............: 1.001
Z-value...................: -1.576594
P-value...................: 0.114889
Hosmer and Lemeshow r^2...: 0.013058
Cox and Snell r^2.........: 0.01627
Nagelkerke's pseudo r^2...: 0.022747
Sample size of AE DB......: 611
Sample size of model......: 600
Missing data %............: 1.800327 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"Z-value", "P-value", "r^2_l", "r^2_cs", "r^2_nagelkerke", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`Z-value` <- as.numeric(GLM.results$`Z-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2_l` <- as.numeric(GLM.results$`r^2_l`)
GLM.results$`r^2_cs` <- as.numeric(GLM.results$`r^2_cs`)
GLM.results$`r^2_nagelkerke` <- as.numeric(GLM.results$`r^2_nagelkerke`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Bin.Uni.Protein.PlaquePhenotypes.RANK.MODEL1.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Bin.Uni.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)In this model we correct for Age, Gender, year of surgery, Hypertension status, Diabetes status, current smoker status, lipid-lowering drugs (LLDs), antiplatelet medication, eGFR (MDRD), BMI, MedHx_CVD (combination of CAD history, stroke history, and peripheral interventions), and stenosis.
Here we use the inverse-rank normalized data - visually
this is more normally distributed.
Analysis of continuous/quantitative plaque traits as a function of plaque HDAC9 expression levels.
GLM.results <- data.frame(matrix(NA, ncol = 15, nrow = 0))
cat("Running linear regression...\n")Running linear regression...for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
for (trait in 1:length(TRAITS.CON.RANK)) {
TRAIT = TRAITS.CON.RANK[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M2) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
### univariate
# fit <- lm(currentDF[,TARGET] ~ currentDF[,TRAIT] + Age + Gender + ORdate_year +
# Hypertension.composite + DiabetesStatus + SmokerStatus +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# MedHx_CVD + stenose,
# data = currentDF)
fit <- lm(currentDF[,TARGET] ~ currentDF[,TRAIT] + Age + Gender + ORdate_epoch +
Hypertension.composite + DiabetesStatus + SmokerStatus +
Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
MedHx_CVD + stenose,
data = currentDF)
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 15, nrow = 0))
GLM.results.TEMP[1,] = GLM.CON(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
Analysis of HDAC9.
- processing MAC_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ 1, data = currentDF)
Coefficients:
(Intercept)
28.44
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-33.50 -17.49 -6.78 7.74 414.80
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.492e+02 2.525e+02 -0.591 0.555
currentDF[, TRAIT] -3.679e-01 1.588e+00 -0.232 0.817
Age 1.993e-01 1.970e-01 1.011 0.312
Gendermale 1.233e+00 3.594e+00 0.343 0.732
ORdate_epoch 1.175e-08 1.880e-08 0.625 0.532
Hypertension.compositeyes 2.040e+00 4.540e+00 0.449 0.653
DiabetesStatusDiabetes -5.176e+00 3.734e+00 -1.386 0.166
SmokerStatusEx-smoker -1.015e-02 3.417e+00 -0.003 0.998
SmokerStatusNever smoked -2.624e+00 4.481e+00 -0.586 0.558
Med.Statin.LLDyes 7.455e-02 3.631e+00 0.021 0.984
Med.all.antiplateletyes 4.478e+00 5.154e+00 0.869 0.385
GFR_MDRD -8.473e-04 2.458e-03 -0.345 0.730
BMI -1.084e-03 5.151e-03 -0.210 0.833
MedHx_CVDNo 7.364e-01 3.140e+00 0.235 0.815
stenose50-70% 1.311e+01 2.467e+01 0.531 0.595
stenose70-90% 7.327e+00 2.400e+01 0.305 0.760
stenose90-99% 1.198e+01 2.402e+01 0.499 0.618
stenose100% (Occlusion) 7.806e+00 2.946e+01 0.265 0.791
Residual standard error: 33.53 on 502 degrees of freedom
Multiple R-squared: 0.01664, Adjusted R-squared: -0.01666
F-statistic: 0.4998 on 17 and 502 DF, p-value: 0.9533
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_rankNorm
Effect size...............: -0.36787
Standard error............: 1.588445
Odds ratio (effect size)..: 0.692
Lower 95% CI..............: 0.031
Upper 95% CI..............: 15.572
T-value...................: -0.231591
P-value...................: 0.81695
R^2.......................: 0.016643
Adjusted r^2..............: -0.016658
Sample size of AE DB......: 611
Sample size of model......: 520
Missing data %............: 14.89362
- processing SMC_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT], data = currentDF)
Coefficients:
(Intercept) currentDF[, TRAIT]
28.147 -5.029
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-40.36 -17.45 -6.66 6.92 413.88
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.197e+01 2.507e+02 -0.207 0.83589
currentDF[, TRAIT] -4.828e+00 1.641e+00 -2.941 0.00342 **
Age 1.512e-01 1.963e-01 0.770 0.44157
Gendermale -4.477e-01 3.602e+00 -0.124 0.90113
ORdate_epoch 4.589e-09 1.865e-08 0.246 0.80573
Hypertension.compositeyes 1.643e+00 4.502e+00 0.365 0.71522
DiabetesStatusDiabetes -4.248e+00 3.711e+00 -1.145 0.25288
SmokerStatusEx-smoker 2.042e-01 3.393e+00 0.060 0.95203
SmokerStatusNever smoked -3.202e+00 4.446e+00 -0.720 0.47177
Med.Statin.LLDyes -3.433e-01 3.610e+00 -0.095 0.92428
Med.all.antiplateletyes 4.117e+00 5.113e+00 0.805 0.42113
GFR_MDRD -1.431e-04 2.448e-03 -0.058 0.95338
BMI -1.149e-03 5.106e-03 -0.225 0.82200
MedHx_CVDNo 1.331e+00 3.121e+00 0.426 0.66996
stenose50-70% 1.345e+01 2.444e+01 0.551 0.58221
stenose70-90% 8.150e+00 2.376e+01 0.343 0.73172
stenose90-99% 1.305e+01 2.380e+01 0.548 0.58373
stenose100% (Occlusion) 9.146e+00 2.922e+01 0.313 0.75445
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33.26 on 501 degrees of freedom
Multiple R-squared: 0.0334, Adjusted R-squared: 0.0005975
F-statistic: 1.018 on 17 and 501 DF, p-value: 0.4362
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' SMC_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: SMC_rankNorm
Effect size...............: -4.82779
Standard error............: 1.641478
Odds ratio (effect size)..: 0.008
Lower 95% CI..............: 0
Upper 95% CI..............: 0.2
T-value...................: -2.941123
P-value...................: 0.003421542
R^2.......................: 0.033396
Adjusted r^2..............: 0.000597
Sample size of AE DB......: 611
Sample size of model......: 519
Missing data %............: 15.05728
- processing MAC_SMC_ratio_rank
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + DiabetesStatus,
data = currentDF)
Coefficients:
(Intercept) currentDF[, TRAIT] DiabetesStatusDiabetes
29.469 -3.705 -5.199
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-37.27 -16.42 -6.72 9.08 410.01
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.416e+02 2.490e+02 -0.970 0.33232
currentDF[, TRAIT] -4.236e+00 1.530e+00 -2.769 0.00583 **
Age 2.064e-01 1.956e-01 1.055 0.29197
Gendermale 6.190e-01 3.570e+00 0.173 0.86240
ORdate_epoch 1.832e-08 1.851e-08 0.990 0.32273
Hypertension.compositeyes 2.418e+00 4.508e+00 0.536 0.59187
DiabetesStatusDiabetes -5.518e+00 3.708e+00 -1.488 0.13737
SmokerStatusEx-smoker -4.992e-01 3.397e+00 -0.147 0.88322
SmokerStatusNever smoked -3.748e+00 4.464e+00 -0.840 0.40154
Med.Statin.LLDyes 2.952e-01 3.607e+00 0.082 0.93481
Med.all.antiplateletyes 5.344e+00 5.129e+00 1.042 0.29797
GFR_MDRD -8.481e-04 2.441e-03 -0.347 0.72837
BMI -1.841e-03 5.118e-03 -0.360 0.71917
MedHx_CVDNo 1.086e+00 3.120e+00 0.348 0.72805
stenose50-70% 1.629e+01 2.449e+01 0.665 0.50626
stenose70-90% 1.116e+01 2.383e+01 0.468 0.63968
stenose90-99% 1.656e+01 2.388e+01 0.694 0.48830
stenose100% (Occlusion) 1.190e+01 2.929e+01 0.406 0.68459
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33.29 on 501 degrees of freedom
Multiple R-squared: 0.03153, Adjusted R-squared: -0.001334
F-statistic: 0.9594 on 17 and 501 DF, p-value: 0.5035
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_SMC_ratio_rank ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_SMC_ratio_rank
Effect size...............: -4.236248
Standard error............: 1.529904
Odds ratio (effect size)..: 0.014
Lower 95% CI..............: 0.001
Upper 95% CI..............: 0.29
T-value...................: -2.768964
P-value...................: 0.00583196
R^2.......................: 0.031528
Adjusted r^2..............: -0.001334
Sample size of AE DB......: 611
Sample size of model......: 519
Missing data %............: 15.05728
- processing VesselDensity_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ DiabetesStatus, data = currentDF)
Coefficients:
(Intercept) DiabetesStatusDiabetes
29.937 -6.025
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-34.56 -18.14 -6.57 8.25 414.27
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.151e+02 2.678e+02 -0.803 0.422
currentDF[, TRAIT] -3.175e-01 1.703e+00 -0.186 0.852
Age 1.910e-01 2.100e-01 0.910 0.363
Gendermale 1.197e+00 3.788e+00 0.316 0.752
ORdate_epoch 1.689e-08 1.992e-08 0.848 0.397
Hypertension.compositeyes 1.661e-01 4.826e+00 0.034 0.973
DiabetesStatusDiabetes -6.119e+00 4.085e+00 -1.498 0.135
SmokerStatusEx-smoker -7.764e-02 3.634e+00 -0.021 0.983
SmokerStatusNever smoked -2.278e+00 4.813e+00 -0.473 0.636
Med.Statin.LLDyes -1.437e+00 3.830e+00 -0.375 0.708
Med.all.antiplateletyes 4.744e+00 5.464e+00 0.868 0.386
GFR_MDRD -1.282e-03 2.610e-03 -0.491 0.624
BMI -2.992e-04 5.456e-03 -0.055 0.956
MedHx_CVDNo 1.580e+00 3.329e+00 0.475 0.635
stenose50-70% 1.230e+01 2.536e+01 0.485 0.628
stenose70-90% 7.525e+00 2.454e+01 0.307 0.759
stenose90-99% 1.221e+01 2.458e+01 0.497 0.620
stenose100% (Occlusion) 7.699e+00 3.017e+01 0.255 0.799
stenose50-99% 2.792e+00 4.273e+01 0.065 0.948
Residual standard error: 34.29 on 468 degrees of freedom
Multiple R-squared: 0.01795, Adjusted R-squared: -0.01982
F-statistic: 0.4752 on 18 and 468 DF, p-value: 0.9679
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' VesselDensity_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: VesselDensity_rankNorm
Effect size...............: -0.317474
Standard error............: 1.702736
Odds ratio (effect size)..: 0.728
Lower 95% CI..............: 0.026
Upper 95% CI..............: 20.489
T-value...................: -0.186449
P-value...................: 0.8521732
R^2.......................: 0.01795
Adjusted r^2..............: -0.019821
Sample size of AE DB......: 611
Sample size of model......: 487
Missing data %............: 20.2946 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"T-value", "P-value", "r^2", "r^2_adj", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`T-value` <- as.numeric(GLM.results$`T-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2` <- as.numeric(GLM.results$`r^2`)
GLM.results$`r^2_adj` <- as.numeric(GLM.results$`r^2_adj`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
library(openxlsx)
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Con.Multi.Protein.PlaquePhenotypes.RANK.MODEL2.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Con.Multi.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)Analysis of binary plaque traits as a function of plaque MCP1 levels.
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
for (trait in 1:length(TRAITS.BIN)) {
TRAIT = TRAITS.BIN[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M2) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_year +
# Hypertension.composite + DiabetesStatus + SmokerStatus +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# MedHx_CVD + stenose,
# data = currentDF, family = binomial(link = "logit"))
fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_epoch +
Hypertension.composite + DiabetesStatus + SmokerStatus +
Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
MedHx_CVD + stenose,
data = currentDF, family = binomial(link = "logit"))
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 16, nrow = 0))
GLM.results.TEMP[1,] = GLM.BIN(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
Analysis of HDAC9.
- processing CalcificationPlaque
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender + ORdate_epoch + Hypertension.composite, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale ORdate_epoch Hypertension.compositeyes
8.113e+01 5.333e-03 -3.216e-01 -6.080e-09 4.737e-01
Degrees of Freedom: 524 Total (i.e. Null); 520 Residual
Null Deviance: 727.5
Residual Deviance: 685.9 AIC: 695.9
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.790 -1.082 0.640 1.074 1.760
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.375e+01 1.561e+01 4.724 2.32e-06 ***
currentDF[, TARGET] 5.116e-03 3.386e-03 1.511 0.1308
Age -6.373e-03 1.225e-02 -0.520 0.6029
Gendermale -3.365e-01 2.242e-01 -1.501 0.1335
ORdate_epoch -5.439e-09 1.161e-09 -4.687 2.77e-06 ***
Hypertension.compositeyes 5.208e-01 2.856e-01 1.824 0.0682 .
DiabetesStatusDiabetes -2.796e-02 2.303e-01 -0.121 0.9034
SmokerStatusEx-smoker -5.353e-02 2.126e-01 -0.252 0.8012
SmokerStatusNever smoked -2.635e-01 2.791e-01 -0.944 0.3450
Med.Statin.LLDyes -1.876e-01 2.252e-01 -0.833 0.4048
Med.all.antiplateletyes -1.568e-01 3.205e-01 -0.489 0.6247
GFR_MDRD -1.717e-04 1.524e-04 -1.127 0.2599
BMI -8.958e-05 3.194e-04 -0.280 0.7792
MedHx_CVDNo 6.958e-03 1.957e-01 0.036 0.9716
stenose50-70% -4.450e-01 1.547e+00 -0.288 0.7737
stenose70-90% -5.208e-01 1.504e+00 -0.346 0.7291
stenose90-99% -1.987e-01 1.506e+00 -0.132 0.8950
stenose100% (Occlusion) 6.063e-01 1.914e+00 0.317 0.7515
stenose50-99% -1.221e+01 5.354e+02 -0.023 0.9818
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 727.48 on 524 degrees of freedom
Residual deviance: 678.49 on 506 degrees of freedom
AIC: 716.49
Number of Fisher Scoring iterations: 12
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' CalcificationPlaque ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: CalcificationPlaque
Effect size...............: 0.005116
Standard error............: 0.003386
Odds ratio (effect size)..: 1.005
Lower 95% CI..............: 0.998
Upper 95% CI..............: 1.012
Z-value...................: 1.51086
P-value...................: 0.1308242
Hosmer and Lemeshow r^2...: 0.06734
Cox and Snell r^2.........: 0.08909
Nagelkerke's pseudo r^2...: 0.118811
Sample size of AE DB......: 611
Sample size of model......: 525
Missing data %............: 14.07529
- processing CollagenPlaque
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ SmokerStatus +
Med.all.antiplatelet, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) SmokerStatusEx-smoker SmokerStatusNever smoked Med.all.antiplateletyes
1.13580 -0.09493 -0.68875 0.56243
Degrees of Freedom: 523 Total (i.e. Null); 520 Residual
Null Deviance: 507.9
Residual Deviance: 499.8 AIC: 507.8
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.1792 0.4865 0.5672 0.6577 1.1449
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.457e+00 1.013e+03 0.003 0.9973
currentDF[, TARGET] -4.143e-03 2.996e-03 -1.383 0.1667
Age 1.631e-02 1.513e-02 1.078 0.2812
Gendermale -6.460e-02 2.790e-01 -0.232 0.8169
ORdate_epoch 8.611e-10 1.427e-09 0.603 0.5462
Hypertension.compositeyes -1.477e-01 3.617e-01 -0.408 0.6830
DiabetesStatusDiabetes 9.145e-02 2.949e-01 0.310 0.7565
SmokerStatusEx-smoker -1.231e-01 2.745e-01 -0.448 0.6539
SmokerStatusNever smoked -7.641e-01 3.241e-01 -2.358 0.0184 *
Med.Statin.LLDyes -1.756e-01 2.857e-01 -0.615 0.5387
Med.all.antiplateletyes 7.610e-01 3.578e-01 2.127 0.0334 *
GFR_MDRD 1.717e-05 1.916e-04 0.090 0.9286
BMI 2.145e-04 3.983e-04 0.539 0.5901
MedHx_CVDNo -7.812e-02 2.431e-01 -0.321 0.7480
stenose50-70% -1.511e+01 1.013e+03 -0.015 0.9881
stenose70-90% -1.431e+01 1.013e+03 -0.014 0.9887
stenose90-99% -1.415e+01 1.013e+03 -0.014 0.9888
stenose100% (Occlusion) 2.109e-01 1.237e+03 0.000 0.9999
stenose50-99% -8.065e-01 1.773e+03 0.000 0.9996
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 507.93 on 523 degrees of freedom
Residual deviance: 488.42 on 505 degrees of freedom
AIC: 526.42
Number of Fisher Scoring iterations: 14
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' CollagenPlaque ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: CollagenPlaque
Effect size...............: -0.004143
Standard error............: 0.002996
Odds ratio (effect size)..: 0.996
Lower 95% CI..............: 0.99
Upper 95% CI..............: 1.002
Z-value...................: -1.382886
P-value...................: 0.1667
Hosmer and Lemeshow r^2...: 0.038422
Cox and Snell r^2.........: 0.036559
Nagelkerke's pseudo r^2...: 0.058903
Sample size of AE DB......: 611
Sample size of model......: 524
Missing data %............: 14.23895
- processing Fat10Perc
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender + ORdate_epoch + BMI, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale ORdate_epoch BMI
9.026e+01 1.994e-02 7.199e-01 -6.702e-09 -5.056e-04
Degrees of Freedom: 524 Total (i.e. Null); 520 Residual
Null Deviance: 580.9
Residual Deviance: 533.2 AIC: 543.2
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.57779 0.07006 0.56602 0.75768 1.36946
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.157e+02 6.173e+02 0.187 0.85134
currentDF[, TARGET] 1.912e-02 6.065e-03 3.152 0.00162 **
Age 1.699e-02 1.443e-02 1.177 0.23903
Gendermale 7.617e-01 2.525e-01 3.016 0.00256 **
ORdate_epoch -7.667e-09 1.379e-09 -5.559 2.72e-08 ***
Hypertension.compositeyes -1.967e-02 3.359e-01 -0.059 0.95331
DiabetesStatusDiabetes -1.488e-01 2.678e-01 -0.556 0.57853
SmokerStatusEx-smoker -1.190e-01 2.513e-01 -0.474 0.63567
SmokerStatusNever smoked -1.620e-02 3.301e-01 -0.049 0.96085
Med.Statin.LLDyes 4.476e-02 2.717e-01 0.165 0.86916
Med.all.antiplateletyes 4.163e-01 3.614e-01 1.152 0.24930
GFR_MDRD 1.224e-04 1.820e-04 0.673 0.50124
BMI -5.080e-04 3.731e-04 -1.362 0.17329
MedHx_CVDNo 3.757e-02 2.298e-01 0.163 0.87015
stenose50-70% -1.266e+01 6.170e+02 -0.021 0.98363
stenose70-90% -1.351e+01 6.170e+02 -0.022 0.98254
stenose90-99% -1.374e+01 6.170e+02 -0.022 0.98224
stenose100% (Occlusion) -1.429e+01 6.170e+02 -0.023 0.98152
stenose50-99% -2.812e+01 1.077e+03 -0.026 0.97917
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 580.93 on 524 degrees of freedom
Residual deviance: 521.97 on 506 degrees of freedom
AIC: 559.97
Number of Fisher Scoring iterations: 13
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' Fat10Perc ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: Fat10Perc
Effect size...............: 0.019116
Standard error............: 0.006065
Odds ratio (effect size)..: 1.019
Lower 95% CI..............: 1.007
Upper 95% CI..............: 1.031
Z-value...................: 3.151719
P-value...................: 0.001623123
Hosmer and Lemeshow r^2...: 0.101483
Cox and Snell r^2.........: 0.106219
Nagelkerke's pseudo r^2...: 0.158703
Sample size of AE DB......: 611
Sample size of model......: 525
Missing data %............: 14.07529
- processing IPH
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender + ORdate_epoch,
family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) Gendermale ORdate_epoch
8.149e+01 4.364e-01 -6.074e-09
Degrees of Freedom: 523 Total (i.e. Null); 521 Residual
Null Deviance: 697.8
Residual Deviance: 662.8 AIC: 668.8
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0537 -1.1708 0.7148 0.9695 1.7294
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 9.619e+01 5.900e+02 0.163 0.8705
currentDF[, TARGET] 7.191e-04 2.932e-03 0.245 0.8063
Age 5.330e-03 1.260e-02 0.423 0.6722
Gendermale 4.544e-01 2.259e-01 2.011 0.0443 *
ORdate_epoch -6.090e-09 1.191e-09 -5.116 3.12e-07 ***
Hypertension.compositeyes -2.525e-01 2.910e-01 -0.868 0.3856
DiabetesStatusDiabetes -1.713e-01 2.357e-01 -0.727 0.4674
SmokerStatusEx-smoker 8.778e-02 2.178e-01 0.403 0.6869
SmokerStatusNever smoked 2.782e-01 2.872e-01 0.969 0.3326
Med.Statin.LLDyes -2.047e-01 2.351e-01 -0.871 0.3839
Med.all.antiplateletyes 2.090e-01 3.266e-01 0.640 0.5222
GFR_MDRD -8.243e-05 1.573e-04 -0.524 0.6004
BMI 2.548e-05 3.279e-04 0.078 0.9381
MedHx_CVDNo -3.298e-01 1.991e-01 -1.656 0.0976 .
stenose50-70% -1.439e+01 5.898e+02 -0.024 0.9805
stenose70-90% -1.443e+01 5.898e+02 -0.024 0.9805
stenose90-99% -1.425e+01 5.898e+02 -0.024 0.9807
stenose100% (Occlusion) -1.512e+01 5.898e+02 -0.026 0.9795
stenose50-99% -2.838e+01 1.062e+03 -0.027 0.9787
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 697.75 on 523 degrees of freedom
Residual deviance: 651.46 on 505 degrees of freedom
AIC: 689.46
Number of Fisher Scoring iterations: 13
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' IPH ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: IPH
Effect size...............: 0.000719
Standard error............: 0.002932
Odds ratio (effect size)..: 1.001
Lower 95% CI..............: 0.995
Upper 95% CI..............: 1.006
Z-value...................: 0.245254
P-value...................: 0.8062599
Hosmer and Lemeshow r^2...: 0.066347
Cox and Snell r^2.........: 0.084556
Nagelkerke's pseudo r^2...: 0.114895
Sample size of AE DB......: 611
Sample size of model......: 524
Missing data %............: 14.23895
- processing MAC_binned
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender + BMI, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) Gendermale BMI
0.2421162 0.5078264 -0.0006331
Degrees of Freedom: 521 Total (i.e. Null); 519 Residual
Null Deviance: 718
Residual Deviance: 708 AIC: 714
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8509 -1.2174 0.8838 1.0686 1.4859
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.769e+01 1.498e+01 1.181 0.2378
currentDF[, TARGET] 3.232e-03 3.085e-03 1.047 0.2949
Age -2.216e-03 1.202e-02 -0.184 0.8537
Gendermale 5.568e-01 2.188e-01 2.544 0.0110 *
ORdate_epoch -1.345e-09 1.115e-09 -1.206 0.2279
Hypertension.compositeyes 8.827e-02 2.748e-01 0.321 0.7480
DiabetesStatusDiabetes -1.574e-01 2.281e-01 -0.690 0.4902
SmokerStatusEx-smoker 9.348e-02 2.087e-01 0.448 0.6542
SmokerStatusNever smoked 1.130e-01 2.731e-01 0.414 0.6791
Med.Statin.LLDyes 2.816e-01 2.202e-01 1.279 0.2010
Med.all.antiplateletyes 2.371e-02 3.137e-01 0.076 0.9397
GFR_MDRD -1.221e-04 1.503e-04 -0.812 0.4167
BMI -7.240e-04 3.146e-04 -2.301 0.0214 *
MedHx_CVDNo -2.032e-02 1.916e-01 -0.106 0.9156
stenose50-70% 1.050e-01 1.476e+00 0.071 0.9433
stenose70-90% 5.171e-01 1.433e+00 0.361 0.7183
stenose90-99% 3.313e-01 1.436e+00 0.231 0.8175
stenose100% (Occlusion) -1.394e-01 1.762e+00 -0.079 0.9370
stenose50-99% -1.298e+01 5.354e+02 -0.024 0.9807
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 718.05 on 521 degrees of freedom
Residual deviance: 697.63 on 503 degrees of freedom
AIC: 735.63
Number of Fisher Scoring iterations: 12
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_binned ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_binned
Effect size...............: 0.003232
Standard error............: 0.003085
Odds ratio (effect size)..: 1.003
Lower 95% CI..............: 0.997
Upper 95% CI..............: 1.009
Z-value...................: 1.047467
P-value...................: 0.2948845
Hosmer and Lemeshow r^2...: 0.02844
Cox and Snell r^2.........: 0.038366
Nagelkerke's pseudo r^2...: 0.051339
Sample size of AE DB......: 611
Sample size of model......: 522
Missing data %............: 14.56629
- processing SMC_binned
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + MedHx_CVD, family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Age MedHx_CVDNo
1.393861 -0.005854 -0.019210 0.454945
Degrees of Freedom: 522 Total (i.e. Null); 519 Residual
Null Deviance: 661.1
Residual Deviance: 648 AIC: 656
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0579 -1.2956 0.7333 0.8952 1.8563
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.058e+01 1.622e+01 -0.652 0.5144
currentDF[, TARGET] -6.555e-03 2.996e-03 -2.188 0.0287 *
Age -2.122e-02 1.312e-02 -1.618 0.1057
Gendermale -3.581e-01 2.428e-01 -1.475 0.1402
ORdate_epoch 8.644e-10 1.210e-09 0.714 0.4752
Hypertension.compositeyes 1.427e-01 2.943e-01 0.485 0.6279
DiabetesStatusDiabetes -2.093e-01 2.396e-01 -0.874 0.3823
SmokerStatusEx-smoker 6.681e-02 2.241e-01 0.298 0.7656
SmokerStatusNever smoked -4.101e-01 2.864e-01 -1.432 0.1521
Med.Statin.LLDyes 5.011e-02 2.354e-01 0.213 0.8314
Med.all.antiplateletyes -1.768e-01 3.374e-01 -0.524 0.6002
GFR_MDRD -4.990e-05 1.600e-04 -0.312 0.7551
BMI 2.333e-04 3.391e-04 0.688 0.4915
MedHx_CVDNo 4.956e-01 2.104e-01 2.356 0.0185 *
stenose50-70% 2.082e-01 1.476e+00 0.141 0.8878
stenose70-90% 5.682e-01 1.432e+00 0.397 0.6916
stenose90-99% 1.037e+00 1.437e+00 0.722 0.4704
stenose100% (Occlusion) -1.405e-01 1.772e+00 -0.079 0.9368
stenose50-99% 1.388e+01 5.354e+02 0.026 0.9793
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 661.08 on 522 degrees of freedom
Residual deviance: 632.23 on 504 degrees of freedom
AIC: 670.23
Number of Fisher Scoring iterations: 12
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' SMC_binned ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: SMC_binned
Effect size...............: -0.006555
Standard error............: 0.002996
Odds ratio (effect size)..: 0.993
Lower 95% CI..............: 0.988
Upper 95% CI..............: 0.999
Z-value...................: -2.187728
P-value...................: 0.02868944
Hosmer and Lemeshow r^2...: 0.04363
Cox and Snell r^2.........: 0.053656
Nagelkerke's pseudo r^2...: 0.074783
Sample size of AE DB......: 611
Sample size of model......: 523
Missing data %............: 14.40262 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"Z-value", "P-value", "r^2_l", "r^2_cs", "r^2_nagelkerke", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`Z-value` <- as.numeric(GLM.results$`Z-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2_l` <- as.numeric(GLM.results$`r^2_l`)
GLM.results$`r^2_cs` <- as.numeric(GLM.results$`r^2_cs`)
GLM.results$`r^2_nagelkerke` <- as.numeric(GLM.results$`r^2_nagelkerke`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Bin.Multi.Protein.PlaquePhenotypes.RANK.MODEL2.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Bin.Multi.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)We will perform a cross-sectional analysis between plaque HDAC9 expression levels and the ‘clinical status’ of the plaque in terms of presence of patients’ symptoms (symptomatic vs. asymptomatic). The symptoms of interest are:
In this model we correct for Age, Gender, and year of surgery.
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
TRAIT = "AsymptSympt"
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M1) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# + Hypertension.composite + DiabetesStatus + SmokerCurrent +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# CAD_history + Stroke_history + Peripheral.interv + stenose
# fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_year,
# data = currentDF, family = binomial(link = "logit"))
fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_epoch,
data = currentDF, family = binomial(link = "logit"))
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 16, nrow = 0))
GLM.results.TEMP[1,] = GLM.BIN(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
Analysis of HDAC9.
- processing AsymptSympt
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender + ORdate_epoch,
family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) Gendermale ORdate_epoch
3.550e+01 3.275e-01 -2.606e-09
Degrees of Freedom: 610 Total (i.e. Null); 608 Residual
Null Deviance: 747.7
Residual Deviance: 738.3 AIC: 744.3
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7976 -1.4056 0.7646 0.8692 1.2170
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.516e+01 1.324e+01 2.656 0.00792 **
currentDF[, TARGET] -3.075e-03 2.614e-03 -1.177 0.23934
Age -3.058e-03 9.998e-03 -0.306 0.75970
Gendermale 3.271e-01 2.023e-01 1.617 0.10595
ORdate_epoch -2.565e-09 9.897e-10 -2.592 0.00954 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 747.66 on 610 degrees of freedom
Residual deviance: 736.86 on 606 degrees of freedom
AIC: 746.86
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' AsymptSympt ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: AsymptSympt
Effect size...............: -0.003075
Standard error............: 0.002614
Odds ratio (effect size)..: 0.997
Lower 95% CI..............: 0.992
Upper 95% CI..............: 1.002
Z-value...................: -1.176627
P-value...................: 0.2393443
Hosmer and Lemeshow r^2...: 0.014449
Cox and Snell r^2.........: 0.017525
Nagelkerke's pseudo r^2...: 0.024829
Sample size of AE DB......: 611
Sample size of model......: 611
Missing data %............: 0 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"Z-value", "P-value", "r^2_l", "r^2_cs", "r^2_nagelkerke", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`Z-value` <- as.numeric(GLM.results$`Z-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2_l` <- as.numeric(GLM.results$`r^2_l`)
GLM.results$`r^2_cs` <- as.numeric(GLM.results$`r^2_cs`)
GLM.results$`r^2_nagelkerke` <- as.numeric(GLM.results$`r^2_nagelkerke`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Bin.Uni.Protein.RANK.Symptoms.MODEL1.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Bin.Uni.Symptoms")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)In this model we correct for Age, Gender, Hypertension status, Diabetes status, current smoker status, lipid-lowering drugs (LLDs), antiplatelet medication, eGFR (MDRD), BMI, MedHx_CVD (combination of CAD history, stroke history, and peripheral interventions), and stenosis..
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
TRAIT = "AsymptSympt"
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M2) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_year +
# Hypertension.composite + DiabetesStatus + SmokerStatus +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# MedHx_CVD + stenose,
# data = currentDF, family = binomial(link = "logit"))
fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_epoch +
Hypertension.composite + DiabetesStatus + SmokerStatus +
Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
MedHx_CVD + stenose,
data = currentDF, family = binomial(link = "logit"))
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 16, nrow = 0))
GLM.results.TEMP[1,] = GLM.BIN(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
Analysis of HDAC9.
- processing AsymptSympt
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender + ORdate_epoch +
SmokerStatus + Med.Statin.LLD, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) Gendermale ORdate_epoch SmokerStatusEx-smoker SmokerStatusNever smoked Med.Statin.LLDyes
4.297e+01 4.081e-01 -3.214e-09 3.289e-01 7.996e-01 3.428e-01
Degrees of Freedom: 531 Total (i.e. Null); 526 Residual
Null Deviance: 665.1
Residual Deviance: 644.5 AIC: 656.5
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0500 -1.2488 0.7118 0.8737 1.4555
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.441e+01 1.518e+01 2.926 0.00343 **
currentDF[, TARGET] -2.286e-03 2.784e-03 -0.821 0.41168
Age -1.548e-02 1.281e-02 -1.209 0.22680
Gendermale 4.700e-01 2.297e-01 2.046 0.04071 *
ORdate_epoch -3.352e-09 1.129e-09 -2.969 0.00298 **
Hypertension.compositeyes 2.935e-01 2.855e-01 1.028 0.30388
DiabetesStatusDiabetes -2.692e-01 2.396e-01 -1.123 0.26129
SmokerStatusEx-smoker 3.799e-01 2.197e-01 1.729 0.08375 .
SmokerStatusNever smoked 8.193e-01 3.076e-01 2.664 0.00773 **
Med.Statin.LLDyes 2.983e-01 2.323e-01 1.284 0.19918
Med.all.antiplateletyes -2.030e-01 3.460e-01 -0.587 0.55743
GFR_MDRD -8.610e-05 1.612e-04 -0.534 0.59326
BMI -1.234e-04 3.333e-04 -0.370 0.71128
MedHx_CVDNo 2.352e-01 2.082e-01 1.130 0.25846
stenose50-70% 8.132e-01 1.480e+00 0.550 0.58256
stenose70-90% 1.334e+00 1.438e+00 0.927 0.35380
stenose90-99% 8.834e-01 1.439e+00 0.614 0.53944
stenose100% (Occlusion) -6.692e-03 1.777e+00 -0.004 0.99700
stenose50-99% 1.517e+01 5.354e+02 0.028 0.97740
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 665.11 on 531 degrees of freedom
Residual deviance: 630.10 on 513 degrees of freedom
AIC: 668.1
Number of Fisher Scoring iterations: 12
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' AsymptSympt ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: AsymptSympt
Effect size...............: -0.002286
Standard error............: 0.002784
Odds ratio (effect size)..: 0.998
Lower 95% CI..............: 0.992
Upper 95% CI..............: 1.003
Z-value...................: -0.820947
P-value...................: 0.4116763
Hosmer and Lemeshow r^2...: 0.052636
Cox and Snell r^2.........: 0.063687
Nagelkerke's pseudo r^2...: 0.089254
Sample size of AE DB......: 611
Sample size of model......: 532
Missing data %............: 12.92962 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"Z-value", "P-value", "r^2_l", "r^2_cs", "r^2_nagelkerke", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`Z-value` <- as.numeric(GLM.results$`Z-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2_l` <- as.numeric(GLM.results$`r^2_l`)
GLM.results$`r^2_cs` <- as.numeric(GLM.results$`r^2_cs`)
GLM.results$`r^2_nagelkerke` <- as.numeric(GLM.results$`r^2_nagelkerke`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Bin.Multi.Protein.RANK.Symptoms.MODEL2.xlsx"),
rowNames = FALSE, colNmes = TRUE, sheetName = "Bin.Multi.Symptoms")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)For the longitudinal analyses of plaque HDAC9 expression levels and secondary cardiovascular events over a three-year follow-up period.
The primary outcome is defined as “a composite of fatal or
non-fatal myocardial infarction, fatal or non-fatal stroke, ruptured
aortic aneurysm, fatal cardiac failure, coronary or peripheral
interventions, leg amputation due to vascular causes, and cardiovascular
death”, i.e. major adverse cardiovascular events (MACE). Variable:
epmajor.3years, these include: - myocardial infarction (MI)
- cerebral infarction (CVA/stroke) - cardiovascular death (exact cause
to be investigated) - cerebral bleeding (CVA/stroke) - fatal myocardial
infarction (MI) - fatal cerebral infarction - fatal cerebral bleeding -
sudden death - fatal heart failure - fatal aneurysm rupture - other
cardiovascular death..
The secondary outcomes will be
epstroke.3years, these include:
epcoronary.3years, these include:
epcvdeath.3years,
these include:
We will use 3-year follow-up, but we will also calculate 30 days and 90 days follow-up ‘time-to-event’ variables. On average there are 365.25 days in a year. We can calculate 30-days and 90-days follow-up time based on the three years follow-up.
cutt.off.30days = (1/365.25) * 30
cutt.off.90days = (1/365.25) * 90
# Fix maximum FU of 30 and 90 days
AEDB.CEA <- AEDB.CEA %>%
mutate(
FU.cutt.off.30days = ifelse(max.followup <= cutt.off.30days, max.followup, cutt.off.30days),
FU.cutt.off.90days = ifelse(max.followup <= cutt.off.90days, max.followup, cutt.off.90days)
)
AEDB.temp <- subset(AEDB.CEA, select = c("STUDY_NUMBER", "Age", "Gender", "Hospital", "Artery_summary",
"max.followup",
"FU.cutt.off.3years",
"FU.cutt.off.30days",
"FU.cutt.off.90days"))
require(labelled)
AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
rm(AEDB.temp)
AERNASE.clin.hdac9 <- merge(AERNASE.clin.hdac9,
subset(AEDB.CEA, select = c("STUDY_NUMBER", "max.followup",
"FU.cutt.off.3years", "FU.cutt.off.30days", "FU.cutt.off.90days",
"ep_major_t_3years", "ep_stroke_t_3years", "ep_coronary_t_3years",
"ep_cvdeath_t_3years",
"epmajor.3years", "epstroke.3years", "epcoronary.3years", "epcvdeath.3years"
)),
by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = TRUE)
AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9, select = c("STUDY_NUMBER", "Age", "Gender", "Hospital", "Artery_summary",
"max.followup",
"FU.cutt.off.3years",
"FU.cutt.off.30days",
"FU.cutt.off.90days"))
require(labelled)
AERNASE.clin.hdac9.temp$Gender <- to_factor(AERNASE.clin.hdac9.temp$Gender)
AERNASE.clin.hdac9.temp$Hospital <- to_factor(AERNASE.clin.hdac9.temp$Hospital)
AERNASE.clin.hdac9.temp$Artery_summary <- to_factor(AERNASE.clin.hdac9.temp$Artery_summary)
DT::datatable(AERNASE.clin.hdac9.temp[1:10,], caption = "Excerpt of the whole AERNASE.clin.hdac9.", rownames = FALSE)
rm(AERNASE.clin.hdac9.temp)Here we will calculate the new 30- and 90-days follow-up of the events and their event-times of interest:
epmajor.3years)epstroke.3years)epcoronary.3years)epcvdeath.3years)avg_days_in_year = 365.25
cutt.off.30days.scaled <- cutt.off.30days * 365.25
cutt.off.90days.scaled <- cutt.off.90days * 365.25
# Event times
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>%
mutate(
ep_major_t_30days = ifelse(ep_major_t_3years * avg_days_in_year <= cutt.off.30days.scaled,
ep_major_t_3years * avg_days_in_year, cutt.off.30days.scaled),
ep_stroke_t_30days = ifelse(ep_stroke_t_3years * avg_days_in_year <= cutt.off.30days.scaled,
ep_stroke_t_3years * avg_days_in_year, cutt.off.30days.scaled),
ep_coronary_t_30days = ifelse(ep_coronary_t_3years * avg_days_in_year <= cutt.off.30days.scaled,
ep_coronary_t_3years * avg_days_in_year, cutt.off.30days.scaled),
ep_cvdeath_t_30days = ifelse(ep_cvdeath_t_3years * avg_days_in_year <= cutt.off.30days.scaled,
ep_cvdeath_t_3years * avg_days_in_year, cutt.off.30days.scaled),
ep_major_t_90days = ifelse(ep_major_t_3years * avg_days_in_year <= cutt.off.90days.scaled,
ep_major_t_3years * avg_days_in_year, cutt.off.90days.scaled),
ep_stroke_t_90days = ifelse(ep_stroke_t_3years * avg_days_in_year <= cutt.off.90days.scaled,
ep_stroke_t_3years * avg_days_in_year, cutt.off.90days.scaled),
ep_coronary_t_90days = ifelse(ep_coronary_t_3years * avg_days_in_year <= cutt.off.90days.scaled,
ep_coronary_t_3years * avg_days_in_year, cutt.off.90days.scaled),
ep_cvdeath_t_90days = ifelse(ep_cvdeath_t_3years * avg_days_in_year <= cutt.off.90days.scaled,
ep_cvdeath_t_3years * avg_days_in_year, cutt.off.90days.scaled)
)
attach(AERNASE.clin.hdac9)
AERNASE.clin.hdac9[,"epmajor.30days"] <- AERNASE.clin.hdac9$epmajor.3years
AERNASE.clin.hdac9$epmajor.30days[epmajor.3years == 1 & ep_major_t_3years > cutt.off.30days] <- 0
AERNASE.clin.hdac9[,"epstroke.30days"] <- AERNASE.clin.hdac9$epstroke.3years
AERNASE.clin.hdac9$epstroke.30days[epstroke.3years == 1 & ep_stroke_t_3years > cutt.off.30days] <- 0
AERNASE.clin.hdac9[,"epcoronary.30days"] <- AERNASE.clin.hdac9$epcoronary.3years
AERNASE.clin.hdac9$epcoronary.30days[epcoronary.3years == 1 & ep_coronary_t_3years > cutt.off.30days] <- 0
AERNASE.clin.hdac9[,"epcvdeath.30days"] <- AERNASE.clin.hdac9$epcvdeath.3years
AERNASE.clin.hdac9$epcvdeath.30days[epcvdeath.3years == 1 & ep_cvdeath_t_3years > cutt.off.30days] <- 0
AERNASE.clin.hdac9[,"epmajor.90days"] <- AERNASE.clin.hdac9$epmajor.3years
AERNASE.clin.hdac9$epmajor.90days[epmajor.3years == 1 & ep_major_t_3years > cutt.off.90days] <- 0
AERNASE.clin.hdac9[,"epstroke.90days"] <- AERNASE.clin.hdac9$epstroke.3years
AERNASE.clin.hdac9$epstroke.90days[epstroke.3years == 1 & ep_stroke_t_3years > cutt.off.90days] <- 0
AERNASE.clin.hdac9[,"epcoronary.90days"] <- AERNASE.clin.hdac9$epcoronary.3years
AERNASE.clin.hdac9$epcoronary.90days[epcoronary.3years == 1 & ep_coronary_t_3years > cutt.off.90days] <- 0
AERNASE.clin.hdac9[,"epcvdeath.90days"] <- AERNASE.clin.hdac9$epcvdeath.3years
AERNASE.clin.hdac9$epcvdeath.90days[epcvdeath.3years == 1 & ep_cvdeath_t_3years > cutt.off.90days] <- 0
detach(AERNASE.clin.hdac9)
AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9, select = c("STUDY_NUMBER", "Age", "Gender", "Hospital", "Artery_summary",
"epmajor.3years", "epstroke.3years", "epcoronary.3years", "epcvdeath.3years",
"epmajor.30days", "epstroke.30days", "epcoronary.30days", "epcvdeath.30days",
"epmajor.90days", "epstroke.90days", "epcoronary.90days", "epcvdeath.90days"))
require(labelled)
AERNASE.clin.hdac9.temp$Gender <- to_factor(AERNASE.clin.hdac9.temp$Gender)
AERNASE.clin.hdac9.temp$Hospital <- to_factor(AERNASE.clin.hdac9.temp$Hospital)
AERNASE.clin.hdac9.temp$Artery_summary <- to_factor(AERNASE.clin.hdac9.temp$Artery_summary)
DT::datatable(AERNASE.clin.hdac9.temp[1:10,], caption = "Excerpt of the whole AERNASE.clin.hdac9.", rownames = FALSE)
rm(AERNASE.clin.hdac9.temp)First we do some sanity checks and inventory the time-to-event and event variables.
# Reference: https://bioconductor.org/packages/devel/bioc/vignettes/MultiAssayExperiment/inst/doc/QuickStartMultiAssay.html
# If you want to suppress warnings and messages when installing/loading packages
# suppressPackageStartupMessages({})
install.packages.auto("survival")
install.packages.auto("survminer")
install.packages.auto("Hmisc")
cat("* Creating function to summarize Cox regression and prepare container for results.")* Creating function to summarize Cox regression and prepare container for results.# Function to get summary statistics from Cox regression model
COX.STAT <- function(coxfit, DATASET, OUTCOME, target_of_interest){
cat("Summarizing Cox regression results for '", target_of_interest ,"' and its association to '",OUTCOME,"' in '",DATASET,"'.\n")
if (nrow(summary(coxfit)$coefficients) == 1) {
output = c(target_of_interest, rep(NA,8))
cat("Model not fitted; probably singular.\n")
}else {
cat("Collecting data.\n\n")
cox.sum <- summary(coxfit)
cox.effectsize = cox.sum$coefficients[1,1]
cox.SE = cox.sum$coefficients[1,3]
cox.HReffect = cox.sum$coefficients[1,2]
cox.CI_low = exp(cox.effectsize - 1.96 * cox.SE)
cox.CI_up = exp(cox.effectsize + 1.96 * cox.SE)
cox.zvalue = cox.sum$coefficients[1,4]
cox.pvalue = cox.sum$coefficients[1,5]
cox.sample_size = cox.sum$n
cox.nevents = cox.sum$nevent
output = c(DATASET, OUTCOME, target_of_interest, cox.effectsize, cox.SE, cox.HReffect, cox.CI_low, cox.CI_up, cox.zvalue, cox.pvalue, cox.sample_size, cox.nevents)
cat("We have collected the following:\n")
cat("Dataset used..............:", DATASET, "\n")
cat("Outcome analyzed..........:", OUTCOME, "\n")
cat("Protein...................:", target_of_interest, "\n")
cat("Effect size...............:", round(cox.effectsize, 6), "\n")
cat("Standard error............:", round(cox.SE, 6), "\n")
cat("Odds ratio (effect size)..:", round(cox.HReffect, 3), "\n")
cat("Lower 95% CI..............:", round(cox.CI_low, 3), "\n")
cat("Upper 95% CI..............:", round(cox.CI_up, 3), "\n")
cat("T-value...................:", round(cox.zvalue, 6), "\n")
cat("P-value...................:", signif(cox.pvalue, 8), "\n")
cat("Sample size in model......:", cox.sample_size, "\n")
cat("Number of events..........:", cox.nevents, "\n")
}
return(output)
print(output)
}
times = c("ep_major_t_3years",
"ep_stroke_t_3years", "ep_coronary_t_3years", "ep_cvdeath_t_3years")
endpoints = c("epmajor.3years",
"epstroke.3years", "epcoronary.3years", "epcvdeath.3years")
cat("* Check the cases per event type - for sanity.")* Check the cases per event type - for sanity.for (events in endpoints){
require(labelled)
print(paste0("Printing the summary of: ",events))
# print(summary(AERNASE.clin.hdac9[,events]))
print(table(AERNASE.clin.hdac9[,events]))
}[1] "Printing the summary of: epmajor.3years"
1 2
529 79
[1] "Printing the summary of: epstroke.3years"
1 2
561 47
[1] "Printing the summary of: epcoronary.3years"
1 2
561 47
[1] "Printing the summary of: epcvdeath.3years"
1 2
581 27 cat("* Check distribution of events over time - for sanity.")* Check distribution of events over time - for sanity.for (eventtimes in times){
print(paste0("Printing the summary of: ",eventtimes))
print(summary(AERNASE.clin.hdac9[,eventtimes]))
}[1] "Printing the summary of: ep_major_t_3years"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 2.912 3.000 2.661 3.000 3.000 3
[1] "Printing the summary of: ep_stroke_t_3years"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 2.929 3.000 2.705 3.000 3.000 3
[1] "Printing the summary of: ep_coronary_t_3years"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 2.923 3.000 2.736 3.000 3.000 3
[1] "Printing the summary of: ep_cvdeath_t_3years"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0137 2.9500 3.0000 2.8238 3.0000 3.0000 3 for (eventtime in times){
print(paste0("Printing the distribution of: ",eventtime))
p <- gghistogram(AERNASE.clin.hdac9, x = eventtime, y = "..count..",
main = eventtime, bins = 15,
xlab = "year", color = uithof_color[16], fill = uithof_color[16], ggtheme = theme_minimal())
print(p)
ggsave(file = paste0(QC_loc, "/",Today,".AERNASE.clin.hdac9.EventDistributionPerYear.",eventtime,".pdf"), plot = last_plot())
}[1] "Printing the distribution of: ep_major_t_3years"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_stroke_t_3years"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_coronary_t_3years"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_cvdeath_t_3years"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).
times30 = c("ep_major_t_30days",
"ep_stroke_t_30days", "ep_coronary_t_30days", "ep_cvdeath_t_30days")
endpoints30 = c("epmajor.30days",
"epstroke.30days", "epcoronary.30days", "epcvdeath.30days")
cat("* Check the cases per event type - for sanity.")* Check the cases per event type - for sanity.for (events in endpoints30){
print(paste0("Printing the summary of: ",events))
# print(summary(AERNASE.clin.hdac9[,events]))
print(table(AERNASE.clin.hdac9[,events]))
}[1] "Printing the summary of: epmajor.30days"
0 2
529 79
[1] "Printing the summary of: epstroke.30days"
0 2
561 47
[1] "Printing the summary of: epcoronary.30days"
0 1 2
559 2 47
[1] "Printing the summary of: epcvdeath.30days"
0 2
581 27 cat("* Check distribution of events over time - for sanity.")* Check distribution of events over time - for sanity.for (eventtimes in times30){
print(paste0("Printing the summary of: ",eventtimes))
print(summary(AERNASE.clin.hdac9[,eventtimes]))
}[1] "Printing the summary of: ep_major_t_30days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 30.00 30.00 29.03 30.00 30.00 3
[1] "Printing the summary of: ep_stroke_t_30days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 30.00 30.00 29.17 30.00 30.00 3
[1] "Printing the summary of: ep_coronary_t_30days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 30.00 30.00 29.67 30.00 30.00 3
[1] "Printing the summary of: ep_cvdeath_t_30days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
5.003 30.000 30.000 29.936 30.000 30.000 3 for (eventtime in times30){
print(paste0("Printing the distribution of: ",eventtime))
p <- gghistogram(AERNASE.clin.hdac9, x = eventtime, y = "..count..",
main = eventtime, bins = 15,
xlab = "days", color = uithof_color[16], fill = uithof_color[16], ggtheme = theme_minimal())
print(p)
ggsave(file = paste0(QC_loc, "/",Today,".AERNASE.clin.hdac9.EventDistributionPer30Days.",eventtime,".pdf"), plot = last_plot())
}[1] "Printing the distribution of: ep_major_t_30days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_stroke_t_30days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_coronary_t_30days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_cvdeath_t_30days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).
times90 = c("ep_major_t_90days",
"ep_stroke_t_90days", "ep_coronary_t_90days", "ep_cvdeath_t_90days")
endpoints90 = c("epmajor.90days",
"epstroke.90days", "epcoronary.90days", "epcvdeath.90days")
cat("* Check the cases per event type - for sanity.")* Check the cases per event type - for sanity.for (events in endpoints90){
print(paste0("Printing the summary of: ",events))
# print(summary(AERNASE.clin.hdac9[,events]))
print(table(AERNASE.clin.hdac9[,events]))
}[1] "Printing the summary of: epmajor.90days"
0 2
529 79
[1] "Printing the summary of: epstroke.90days"
0 1 2
560 1 47
[1] "Printing the summary of: epcoronary.90days"
0 1 2
559 2 47
[1] "Printing the summary of: epcvdeath.90days"
0 2
581 27 cat("* Check distribution of events over time - for sanity.")* Check distribution of events over time - for sanity.for (eventtimes in times90){
print(paste0("Printing the summary of: ",eventtimes))
print(summary(AERNASE.clin.hdac9[,eventtimes]))
}[1] "Printing the summary of: ep_major_t_90days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 90.00 90.00 86.73 90.00 90.00 3
[1] "Printing the summary of: ep_stroke_t_90days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 90.00 90.00 87.19 90.00 90.00 3
[1] "Printing the summary of: ep_coronary_t_90days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 90.00 90.00 88.73 90.00 90.00 3
[1] "Printing the summary of: ep_cvdeath_t_90days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
5.003 90.000 90.000 89.668 90.000 90.000 3 for (eventtime in times90){
print(paste0("Printing the distribution of: ",eventtime))
p <- gghistogram(AERNASE.clin.hdac9, x = eventtime, y = "..count..",
main = eventtime, bins = 15,
xlab = "days", color = uithof_color[16], fill = uithof_color[16], ggtheme = theme_minimal())
print(p)
ggsave(file = paste0(QC_loc, "/",Today,".AERNASE.clin.hdac9.EventDistributionPer90Days.",eventtime,".pdf"), plot = last_plot())
}[1] "Printing the distribution of: ep_major_t_90days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_stroke_t_90days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_coronary_t_90days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_cvdeath_t_90days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).
Let’s perform the actual Cox-regressions. We will apply a couple of models:
# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times)){
eptime = times[i]
ep = endpoints[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [years]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 1 (Simple model)
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age+Gender + ORdate_year, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age+Gender + ORdate_year, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
# ylim = c(0.2, 1), xlim = c(0,3), col = c("#595A5C", "#DB003F", "#1290D9"),
ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [years]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL1.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}* Analyzing the effect of plaque target-of-interest on [epmajor.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
314 297
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 14,15,16,17 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 13,14,15,16 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year, data = TEMP.DF)
n= 608, number of events= 79
(3 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -3.570e-01 6.998e-01 2.348e-01 -1.520 0.1285
Age 3.071e-02 1.031e+00 1.378e-02 2.229 0.0258 *
Gendermale 7.655e-01 2.150e+00 3.277e-01 2.336 0.0195 *
ORdate_year2003 -1.131e+00 3.228e-01 5.869e-01 -1.927 0.0540 .
ORdate_year2004 -7.715e-01 4.623e-01 5.189e-01 -1.487 0.1371
ORdate_year2005 -4.270e-01 6.524e-01 4.701e-01 -0.908 0.3637
ORdate_year2006 -3.803e-01 6.836e-01 4.692e-01 -0.811 0.4176
ORdate_year2007 -1.480e+00 2.275e-01 6.281e-01 -2.357 0.0184 *
ORdate_year2008 -2.539e-01 7.758e-01 4.763e-01 -0.533 0.5940
ORdate_year2009 -9.296e-01 3.947e-01 5.392e-01 -1.724 0.0847 .
ORdate_year2010 -5.937e-01 5.523e-01 5.869e-01 -1.012 0.3117
ORdate_year2011 -1.509e+00 2.211e-01 8.052e-01 -1.874 0.0609 .
ORdate_year2012 -8.457e-01 4.293e-01 6.949e-01 -1.217 0.2236
ORdate_year2013 -1.715e+01 3.569e-08 3.792e+03 -0.005 0.9964
ORdate_year2014 -1.681e+01 4.994e-08 4.843e+03 -0.003 0.9972
ORdate_year2015 -1.612e+01 9.973e-08 4.743e+03 -0.003 0.9973
ORdate_year2016 -1.630e+01 8.375e-08 8.557e+03 -0.002 0.9985
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
ORdate_year2020 NA NA 0.000e+00 NA NA
ORdate_year2021 NA NA 0.000e+00 NA NA
ORdate_year2022 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 6.998e-01 1.429e+00 0.44164 1.1088
Age 1.031e+00 9.698e-01 1.00372 1.0594
Gendermale 2.150e+00 4.651e-01 1.13108 4.0874
ORdate_year2003 3.228e-01 3.098e+00 0.10219 1.0197
ORdate_year2004 4.623e-01 2.163e+00 0.16720 1.2783
ORdate_year2005 6.524e-01 1.533e+00 0.25967 1.6394
ORdate_year2006 6.836e-01 1.463e+00 0.27255 1.7147
ORdate_year2007 2.275e-01 4.395e+00 0.06643 0.7793
ORdate_year2008 7.758e-01 1.289e+00 0.30498 1.9732
ORdate_year2009 3.947e-01 2.533e+00 0.13720 1.1357
ORdate_year2010 5.523e-01 1.811e+00 0.17480 1.7448
ORdate_year2011 2.211e-01 4.524e+00 0.04561 1.0713
ORdate_year2012 4.293e-01 2.330e+00 0.10997 1.6757
ORdate_year2013 3.569e-08 2.802e+07 0.00000 Inf
ORdate_year2014 4.994e-08 2.002e+07 0.00000 Inf
ORdate_year2015 9.973e-08 1.003e+07 0.00000 Inf
ORdate_year2016 8.375e-08 1.194e+07 0.00000 Inf
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
ORdate_year2020 NA NA NA NA
ORdate_year2021 NA NA NA NA
ORdate_year2022 NA NA NA NA
Concordance= 0.659 (se = 0.029 )
Likelihood ratio test= 29.34 on 17 df, p=0.03
Wald test = 23.89 on 17 df, p=0.1
Score (logrank) test = 26.91 on 17 df, p=0.06
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epmajor.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epmajor.3years
Protein...................: HDAC9
Effect size...............: -0.35699
Standard error............: 0.234834
Odds ratio (effect size)..: 0.7
Lower 95% CI..............: 0.442
Upper 95% CI..............: 1.109
T-value...................: -1.520178
P-value...................: 0.1284662
Sample size in model......: 608
Number of events..........: 79
* Analyzing the effect of plaque target-of-interest on [epstroke.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
314 297
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 14,15,16,17 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 13,14,15,16 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year, data = TEMP.DF)
n= 608, number of events= 47
(3 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -3.767e-01 6.862e-01 3.043e-01 -1.238 0.2157
Age 2.715e-02 1.028e+00 1.774e-02 1.530 0.1259
Gendermale 3.733e-01 1.453e+00 3.746e-01 0.997 0.3190
ORdate_year2003 -9.921e-01 3.708e-01 7.657e-01 -1.296 0.1951
ORdate_year2004 -8.059e-01 4.467e-01 7.090e-01 -1.137 0.2557
ORdate_year2005 -3.037e-01 7.381e-01 6.135e-01 -0.495 0.6206
ORdate_year2006 -7.952e-03 9.921e-01 5.921e-01 -0.013 0.9893
ORdate_year2007 -1.522e+00 2.182e-01 8.675e-01 -1.755 0.0793 .
ORdate_year2008 -1.811e-01 8.344e-01 6.274e-01 -0.289 0.7729
ORdate_year2009 -1.517e+00 2.194e-01 8.714e-01 -1.741 0.0817 .
ORdate_year2010 -1.578e-01 8.540e-01 7.092e-01 -0.222 0.8239
ORdate_year2011 -1.554e+00 2.113e-01 1.122e+00 -1.385 0.1660
ORdate_year2012 -5.545e-01 5.744e-01 8.714e-01 -0.636 0.5246
ORdate_year2013 -1.699e+01 4.164e-08 5.061e+03 -0.003 0.9973
ORdate_year2014 -1.666e+01 5.831e-08 6.455e+03 -0.003 0.9979
ORdate_year2015 -1.622e+01 9.003e-08 6.501e+03 -0.002 0.9980
ORdate_year2016 -1.654e+01 6.526e-08 1.140e+04 -0.001 0.9988
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
ORdate_year2020 NA NA 0.000e+00 NA NA
ORdate_year2021 NA NA 0.000e+00 NA NA
ORdate_year2022 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 6.862e-01 1.457e+00 0.37795 1.246
Age 1.028e+00 9.732e-01 0.99241 1.064
Gendermale 1.453e+00 6.884e-01 0.69706 3.027
ORdate_year2003 3.708e-01 2.697e+00 0.08268 1.663
ORdate_year2004 4.467e-01 2.239e+00 0.11131 1.793
ORdate_year2005 7.381e-01 1.355e+00 0.22175 2.457
ORdate_year2006 9.921e-01 1.008e+00 0.31087 3.166
ORdate_year2007 2.182e-01 4.582e+00 0.03985 1.195
ORdate_year2008 8.344e-01 1.198e+00 0.24397 2.854
ORdate_year2009 2.194e-01 4.558e+00 0.03976 1.210
ORdate_year2010 8.540e-01 1.171e+00 0.21271 3.429
ORdate_year2011 2.113e-01 4.732e+00 0.02344 1.906
ORdate_year2012 5.744e-01 1.741e+00 0.10409 3.169
ORdate_year2013 4.164e-08 2.402e+07 0.00000 Inf
ORdate_year2014 5.831e-08 1.715e+07 0.00000 Inf
ORdate_year2015 9.003e-08 1.111e+07 0.00000 Inf
ORdate_year2016 6.526e-08 1.532e+07 0.00000 Inf
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
ORdate_year2020 NA NA NA NA
ORdate_year2021 NA NA NA NA
ORdate_year2022 NA NA NA NA
Concordance= 0.678 (se = 0.035 )
Likelihood ratio test= 19.58 on 17 df, p=0.3
Wald test = 15.36 on 17 df, p=0.6
Score (logrank) test = 18.04 on 17 df, p=0.4
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epstroke.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epstroke.3years
Protein...................: HDAC9
Effect size...............: -0.376659
Standard error............: 0.304261
Odds ratio (effect size)..: 0.686
Lower 95% CI..............: 0.378
Upper 95% CI..............: 1.246
T-value...................: -1.237946
P-value...................: 0.215736
Sample size in model......: 608
Number of events..........: 47
* Analyzing the effect of plaque target-of-interest on [epcoronary.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
314 297
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 14,15,16,17 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 13,14,15,16 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year, data = TEMP.DF)
n= 608, number of events= 47
(3 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -3.436e-01 7.092e-01 3.018e-01 -1.138 0.2549
Age 5.263e-03 1.005e+00 1.746e-02 0.301 0.7631
Gendermale 1.087e+00 2.966e+00 4.758e-01 2.285 0.0223 *
ORdate_year2003 -6.180e-01 5.390e-01 6.069e-01 -1.018 0.3086
ORdate_year2004 -1.436e+00 2.379e-01 7.325e-01 -1.960 0.0500 *
ORdate_year2005 -6.738e-01 5.097e-01 5.879e-01 -1.146 0.2517
ORdate_year2006 -3.341e-01 7.160e-01 5.482e-01 -0.609 0.5423
ORdate_year2007 -1.082e+00 3.389e-01 6.726e-01 -1.609 0.1076
ORdate_year2008 -7.672e-01 4.643e-01 6.336e-01 -1.211 0.2260
ORdate_year2009 -1.078e+00 3.401e-01 6.741e-01 -1.600 0.1096
ORdate_year2010 -1.868e+00 1.545e-01 1.096e+00 -1.704 0.0884 .
ORdate_year2011 -1.777e+00 1.691e-01 1.098e+00 -1.619 0.1055
ORdate_year2012 -1.586e+00 2.047e-01 1.100e+00 -1.442 0.1494
ORdate_year2013 -1.739e+01 2.809e-08 5.091e+03 -0.003 0.9973
ORdate_year2014 -1.714e+01 3.587e-08 6.484e+03 -0.003 0.9979
ORdate_year2015 -1.643e+01 7.331e-08 6.318e+03 -0.003 0.9979
ORdate_year2016 -1.629e+01 8.444e-08 1.132e+04 -0.001 0.9989
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
ORdate_year2020 NA NA 0.000e+00 NA NA
ORdate_year2021 NA NA 0.000e+00 NA NA
ORdate_year2022 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 7.092e-01 1.410e+00 0.39254 1.2814
Age 1.005e+00 9.948e-01 0.97146 1.0403
Gendermale 2.966e+00 3.371e-01 1.16735 7.5375
ORdate_year2003 5.390e-01 1.855e+00 0.16405 1.7710
ORdate_year2004 2.379e-01 4.204e+00 0.05661 0.9997
ORdate_year2005 5.097e-01 1.962e+00 0.16104 1.6135
ORdate_year2006 7.160e-01 1.397e+00 0.24450 2.0967
ORdate_year2007 3.389e-01 2.951e+00 0.09069 1.2663
ORdate_year2008 4.643e-01 2.154e+00 0.13413 1.6075
ORdate_year2009 3.401e-01 2.940e+00 0.09076 1.2748
ORdate_year2010 1.545e-01 6.474e+00 0.01802 1.3241
ORdate_year2011 1.691e-01 5.914e+00 0.01966 1.4544
ORdate_year2012 2.047e-01 4.884e+00 0.02371 1.7684
ORdate_year2013 2.809e-08 3.560e+07 0.00000 Inf
ORdate_year2014 3.587e-08 2.787e+07 0.00000 Inf
ORdate_year2015 7.331e-08 1.364e+07 0.00000 Inf
ORdate_year2016 8.444e-08 1.184e+07 0.00000 Inf
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
ORdate_year2020 NA NA NA NA
ORdate_year2021 NA NA NA NA
ORdate_year2022 NA NA NA NA
Concordance= 0.668 (se = 0.037 )
Likelihood ratio test= 19.7 on 17 df, p=0.3
Wald test = 15.53 on 17 df, p=0.6
Score (logrank) test = 18.03 on 17 df, p=0.4
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epcoronary.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epcoronary.3years
Protein...................: HDAC9
Effect size...............: -0.343584
Standard error............: 0.301808
Odds ratio (effect size)..: 0.709
Lower 95% CI..............: 0.393
Upper 95% CI..............: 1.281
T-value...................: -1.138421
P-value...................: 0.2549448
Sample size in model......: 608
Number of events..........: 47
* Analyzing the effect of plaque target-of-interest on [epcvdeath.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
314 297
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 12,13,14,15,16,17 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 11,12,13,14,15,16 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year, data = TEMP.DF)
n= 608, number of events= 27
(3 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -8.607e-01 4.229e-01 4.346e-01 -1.980 0.04768 *
Age 7.113e-02 1.074e+00 2.628e-02 2.707 0.00679 **
Gendermale 6.115e-01 1.843e+00 5.456e-01 1.121 0.26237
ORdate_year2003 -1.596e+00 2.027e-01 8.395e-01 -1.901 0.05732 .
ORdate_year2004 -1.704e+00 1.819e-01 8.376e-01 -2.035 0.04187 *
ORdate_year2005 -9.918e-01 3.709e-01 6.332e-01 -1.566 0.11730
ORdate_year2006 -1.827e+00 1.609e-01 8.390e-01 -2.178 0.02942 *
ORdate_year2007 -1.713e+00 1.804e-01 8.386e-01 -2.042 0.04111 *
ORdate_year2008 -1.245e+00 2.878e-01 7.311e-01 -1.703 0.08850 .
ORdate_year2009 -9.372e-01 3.917e-01 6.856e-01 -1.367 0.17163
ORdate_year2010 -1.059e+00 3.469e-01 8.396e-01 -1.261 0.20731
ORdate_year2011 -1.917e+01 4.712e-09 5.489e+03 -0.003 0.99721
ORdate_year2012 -1.930e+01 4.137e-09 6.718e+03 -0.003 0.99771
ORdate_year2013 -1.960e+01 3.087e-09 1.480e+04 -0.001 0.99894
ORdate_year2014 -1.877e+01 7.060e-09 1.853e+04 -0.001 0.99919
ORdate_year2015 -1.790e+01 1.677e-08 1.934e+04 -0.001 0.99926
ORdate_year2016 -1.875e+01 7.211e-09 3.554e+04 -0.001 0.99958
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
ORdate_year2020 NA NA 0.000e+00 NA NA
ORdate_year2021 NA NA 0.000e+00 NA NA
ORdate_year2022 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 4.229e-01 2.365e+00 0.18041 0.9913
Age 1.074e+00 9.313e-01 1.01982 1.1305
Gendermale 1.843e+00 5.425e-01 0.63265 5.3700
ORdate_year2003 2.027e-01 4.932e+00 0.03912 1.0509
ORdate_year2004 1.819e-01 5.497e+00 0.03523 0.9393
ORdate_year2005 3.709e-01 2.696e+00 0.10722 1.2832
ORdate_year2006 1.609e-01 6.216e+00 0.03107 0.8329
ORdate_year2007 1.804e-01 5.544e+00 0.03486 0.9332
ORdate_year2008 2.878e-01 3.474e+00 0.06867 1.2063
ORdate_year2009 3.917e-01 2.553e+00 0.10219 1.5016
ORdate_year2010 3.469e-01 2.883e+00 0.06691 1.7984
ORdate_year2011 4.712e-09 2.122e+08 0.00000 Inf
ORdate_year2012 4.137e-09 2.417e+08 0.00000 Inf
ORdate_year2013 3.087e-09 3.239e+08 0.00000 Inf
ORdate_year2014 7.060e-09 1.417e+08 0.00000 Inf
ORdate_year2015 1.677e-08 5.964e+07 0.00000 Inf
ORdate_year2016 7.211e-09 1.387e+08 0.00000 Inf
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
ORdate_year2020 NA NA NA NA
ORdate_year2021 NA NA NA NA
ORdate_year2022 NA NA NA NA
Concordance= 0.792 (se = 0.037 )
Likelihood ratio test= 29.29 on 17 df, p=0.03
Wald test = 15.29 on 17 df, p=0.6
Score (logrank) test = 29.48 on 17 df, p=0.03
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epcvdeath.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epcvdeath.3years
Protein...................: HDAC9
Effect size...............: -0.860664
Standard error............: 0.434641
Odds ratio (effect size)..: 0.423
Lower 95% CI..............: 0.18
Upper 95% CI..............: 0.991
T-value...................: -1.980172
P-value...................: 0.04768416
Sample size in model......: 608
Number of events..........: 27 
cat("- Edit the column names...\n")- Edit the column names...colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")- Correct the variable types...COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
cat("- Writing results to Excel-file...\n")- Writing results to Excel-file...head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL1.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")- Removing intermediate files...rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times)){
eptime = times[i]
ep = endpoints[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
#as.integer(TEMP.DF[,ep] == "Excluded")
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [years]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 2 adjusted for age, sex, hypertension, diabetes, smoking, LDL-C levels, lipid-lowering drugs, antiplatelet drugs, eGFR, BMI, history of CVD, level of stenosis
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
# ylim = c(0.2, 1), xlim = c(0,3), col = c("#595A5C", "#DB003F", "#1290D9"),
ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [years]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL2.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}* Analyzing the effect of plaque target-of-interest on [epmajor.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
314 297
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 14,15,16,17,33,34,35 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 13,14,15,16,32,33,34 ; coefficient may be infinite.
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = TEMP.DF)
n= 530, number of events= 67
(81 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -1.318e-01 8.765e-01 2.604e-01 -0.506 0.61280
Age 4.473e-02 1.046e+00 1.822e-02 2.454 0.01411 *
Gendermale 7.081e-01 2.030e+00 3.675e-01 1.927 0.05401 .
ORdate_year2003 -1.039e+00 3.540e-01 6.559e-01 -1.583 0.11334
ORdate_year2004 -5.447e-01 5.800e-01 5.701e-01 -0.955 0.33937
ORdate_year2005 -6.583e-02 9.363e-01 5.329e-01 -0.124 0.90168
ORdate_year2006 6.261e-03 1.006e+00 5.220e-01 0.012 0.99043
ORdate_year2007 -1.106e+00 3.309e-01 6.744e-01 -1.640 0.10107
ORdate_year2008 7.153e-02 1.074e+00 5.518e-01 0.130 0.89686
ORdate_year2009 -9.449e-01 3.887e-01 6.423e-01 -1.471 0.14124
ORdate_year2010 -2.599e-01 7.711e-01 6.334e-01 -0.410 0.68158
ORdate_year2011 -1.020e+00 3.605e-01 8.449e-01 -1.208 0.22715
ORdate_year2012 -2.814e-01 7.548e-01 8.582e-01 -0.328 0.74302
ORdate_year2013 -1.658e+01 6.296e-08 4.559e+03 -0.004 0.99710
ORdate_year2014 -1.612e+01 9.973e-08 1.005e+04 -0.002 0.99872
ORdate_year2015 -1.616e+01 9.562e-08 4.869e+03 -0.003 0.99735
ORdate_year2016 -1.527e+01 2.330e-07 1.005e+04 -0.002 0.99879
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
ORdate_year2020 NA NA 0.000e+00 NA NA
ORdate_year2021 NA NA 0.000e+00 NA NA
ORdate_year2022 NA NA 0.000e+00 NA NA
Hypertension.compositeyes 7.024e-01 2.019e+00 5.423e-01 1.295 0.19529
DiabetesStatusDiabetes 7.616e-01 2.142e+00 2.689e-01 2.832 0.00463 **
SmokerStatusEx-smoker -5.310e-01 5.880e-01 2.720e-01 -1.952 0.05093 .
SmokerStatusNever smoked -8.937e-01 4.091e-01 4.383e-01 -2.039 0.04143 *
Med.Statin.LLDyes -2.653e-01 7.670e-01 2.906e-01 -0.913 0.36141
Med.all.antiplateletyes -2.405e-01 7.863e-01 4.075e-01 -0.590 0.55513
GFR_MDRD -1.692e-04 9.998e-01 1.982e-04 -0.853 0.39341
BMI 7.750e-04 1.001e+00 4.414e-04 1.756 0.07916 .
MedHx_CVDNo -6.775e-01 5.079e-01 3.037e-01 -2.231 0.02569 *
stenose50-70% 1.549e+01 5.346e+06 5.017e+03 0.003 0.99754
stenose70-90% 1.611e+01 9.964e+06 5.017e+03 0.003 0.99744
stenose90-99% 1.618e+01 1.060e+07 5.017e+03 0.003 0.99743
stenose100% (Occlusion) 5.304e-01 1.700e+00 6.605e+03 0.000 0.99994
stenoseNA NA NA 0.000e+00 NA NA
stenose50-99% -6.599e-01 5.169e-01 1.123e+04 0.000 0.99995
stenose70-99% NA NA 0.000e+00 NA NA
stenose99 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 8.765e-01 1.141e+00 0.52619 1.4602
Age 1.046e+00 9.563e-01 1.00905 1.0838
Gendermale 2.030e+00 4.926e-01 0.98786 4.1723
ORdate_year2003 3.540e-01 2.825e+00 0.09786 1.2802
ORdate_year2004 5.800e-01 1.724e+00 0.18974 1.7730
ORdate_year2005 9.363e-01 1.068e+00 0.32949 2.6606
ORdate_year2006 1.006e+00 9.938e-01 0.36174 2.7992
ORdate_year2007 3.309e-01 3.022e+00 0.08823 1.2411
ORdate_year2008 1.074e+00 9.310e-01 0.36421 3.1679
ORdate_year2009 3.887e-01 2.573e+00 0.11039 1.3688
ORdate_year2010 7.711e-01 1.297e+00 0.22282 2.6687
ORdate_year2011 3.605e-01 2.774e+00 0.06881 1.8880
ORdate_year2012 7.548e-01 1.325e+00 0.14038 4.0579
ORdate_year2013 6.296e-08 1.588e+07 0.00000 Inf
ORdate_year2014 9.973e-08 1.003e+07 0.00000 Inf
ORdate_year2015 9.562e-08 1.046e+07 0.00000 Inf
ORdate_year2016 2.329e-07 4.293e+06 0.00000 Inf
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
ORdate_year2020 NA NA NA NA
ORdate_year2021 NA NA NA NA
ORdate_year2022 NA NA NA NA
Hypertension.compositeyes 2.019e+00 4.954e-01 0.69727 5.8437
DiabetesStatusDiabetes 2.142e+00 4.669e-01 1.26426 3.6279
SmokerStatusEx-smoker 5.880e-01 1.701e+00 0.34503 1.0022
SmokerStatusNever smoked 4.091e-01 2.444e+00 0.17331 0.9659
Med.Statin.LLDyes 7.670e-01 1.304e+00 0.43391 1.3558
Med.all.antiplateletyes 7.863e-01 1.272e+00 0.35377 1.7475
GFR_MDRD 9.998e-01 1.000e+00 0.99944 1.0002
BMI 1.001e+00 9.992e-01 0.99991 1.0016
MedHx_CVDNo 5.079e-01 1.969e+00 0.28006 0.9210
stenose50-70% 5.346e+06 1.871e-07 0.00000 Inf
stenose70-90% 9.964e+06 1.004e-07 0.00000 Inf
stenose90-99% 1.060e+07 9.435e-08 0.00000 Inf
stenose100% (Occlusion) 1.700e+00 5.883e-01 0.00000 Inf
stenoseNA NA NA NA NA
stenose50-99% 5.169e-01 1.935e+00 0.00000 Inf
stenose70-99% NA NA NA NA
stenose99 NA NA NA NA
Concordance= 0.763 (se = 0.026 )
Likelihood ratio test= 59 on 31 df, p=0.002
Wald test = 44.52 on 31 df, p=0.05
Score (logrank) test = 57.21 on 31 df, p=0.003
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epmajor.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epmajor.3years
Protein...................: HDAC9
Effect size...............: -0.131769
Standard error............: 0.260373
Odds ratio (effect size)..: 0.877
Lower 95% CI..............: 0.526
Upper 95% CI..............: 1.46
T-value...................: -0.506077
P-value...................: 0.6128024
Sample size in model......: 530
Number of events..........: 67
* Analyzing the effect of plaque target-of-interest on [epstroke.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
314 297
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 14,15,16,17,33,34,35 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 13,14,15,16,32,33,34 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = TEMP.DF)
n= 530, number of events= 37
(81 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -1.703e-01 8.434e-01 3.510e-01 -0.485 0.6276
Age 5.147e-02 1.053e+00 2.370e-02 2.172 0.0298 *
Gendermale 2.527e-01 1.287e+00 4.348e-01 0.581 0.5612
ORdate_year2003 -1.044e+00 3.522e-01 9.289e-01 -1.124 0.2612
ORdate_year2004 -6.143e-01 5.410e-01 8.312e-01 -0.739 0.4599
ORdate_year2005 -1.172e-01 8.894e-01 7.372e-01 -0.159 0.8736
ORdate_year2006 4.584e-01 1.582e+00 6.917e-01 0.663 0.5075
ORdate_year2007 -1.053e+00 3.487e-01 9.485e-01 -1.111 0.2667
ORdate_year2008 7.231e-02 1.075e+00 7.590e-01 0.095 0.9241
ORdate_year2009 -1.800e+00 1.654e-01 1.192e+00 -1.510 0.1311
ORdate_year2010 1.386e-01 1.149e+00 8.076e-01 0.172 0.8638
ORdate_year2011 -9.355e-01 3.924e-01 1.181e+00 -0.792 0.4283
ORdate_year2012 -1.095e-01 8.963e-01 1.209e+00 -0.091 0.9278
ORdate_year2013 -1.640e+01 7.545e-08 6.617e+03 -0.002 0.9980
ORdate_year2014 -1.597e+01 1.158e-07 1.419e+04 -0.001 0.9991
ORdate_year2015 -1.627e+01 8.581e-08 7.550e+03 -0.002 0.9983
ORdate_year2016 -1.526e+01 2.347e-07 1.419e+04 -0.001 0.9991
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
ORdate_year2020 NA NA 0.000e+00 NA NA
ORdate_year2021 NA NA 0.000e+00 NA NA
ORdate_year2022 NA NA 0.000e+00 NA NA
Hypertension.compositeyes 3.872e-01 1.473e+00 6.468e-01 0.599 0.5494
DiabetesStatusDiabetes 5.080e-01 1.662e+00 3.825e-01 1.328 0.1842
SmokerStatusEx-smoker -5.848e-01 5.572e-01 3.580e-01 -1.634 0.1023
SmokerStatusNever smoked -1.694e+00 1.837e-01 7.594e-01 -2.231 0.0257 *
Med.Statin.LLDyes -2.400e-01 7.866e-01 3.863e-01 -0.621 0.5344
Med.all.antiplateletyes -1.004e-01 9.045e-01 5.734e-01 -0.175 0.8610
GFR_MDRD 2.910e-04 1.000e+00 2.682e-04 1.085 0.2779
BMI 1.079e-03 1.001e+00 5.996e-04 1.799 0.0721 .
MedHx_CVDNo -5.917e-01 5.534e-01 3.977e-01 -1.488 0.1368
stenose50-70% 1.563e+01 6.126e+06 6.652e+03 0.002 0.9981
stenose70-90% 1.568e+01 6.437e+06 6.652e+03 0.002 0.9981
stenose90-99% 1.574e+01 6.875e+06 6.652e+03 0.002 0.9981
stenose100% (Occlusion) -1.606e-01 8.516e-01 8.497e+03 0.000 1.0000
stenoseNA NA NA 0.000e+00 NA NA
stenose50-99% -1.751e+00 1.737e-01 1.567e+04 0.000 0.9999
stenose70-99% NA NA 0.000e+00 NA NA
stenose99 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 8.434e-01 1.186e+00 0.42391 1.6782
Age 1.053e+00 9.498e-01 1.00505 1.1029
Gendermale 1.287e+00 7.767e-01 0.54902 3.0190
ORdate_year2003 3.522e-01 2.840e+00 0.05703 2.1747
ORdate_year2004 5.410e-01 1.848e+00 0.10609 2.7591
ORdate_year2005 8.894e-01 1.124e+00 0.20968 3.7722
ORdate_year2006 1.582e+00 6.323e-01 0.40766 6.1356
ORdate_year2007 3.487e-01 2.868e+00 0.05434 2.2380
ORdate_year2008 1.075e+00 9.302e-01 0.24287 4.7580
ORdate_year2009 1.654e-01 6.047e+00 0.01599 1.7100
ORdate_year2010 1.149e+00 8.706e-01 0.23590 5.5927
ORdate_year2011 3.924e-01 2.548e+00 0.03876 3.9720
ORdate_year2012 8.963e-01 1.116e+00 0.08386 9.5788
ORdate_year2013 7.545e-08 1.325e+07 0.00000 Inf
ORdate_year2014 1.158e-07 8.636e+06 0.00000 Inf
ORdate_year2015 8.581e-08 1.165e+07 0.00000 Inf
ORdate_year2016 2.347e-07 4.261e+06 0.00000 Inf
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
ORdate_year2020 NA NA NA NA
ORdate_year2021 NA NA NA NA
ORdate_year2022 NA NA NA NA
Hypertension.compositeyes 1.473e+00 6.789e-01 0.41460 5.2326
DiabetesStatusDiabetes 1.662e+00 6.017e-01 0.78520 3.5174
SmokerStatusEx-smoker 5.572e-01 1.795e+00 0.27625 1.1239
SmokerStatusNever smoked 1.837e-01 5.443e+00 0.04147 0.8139
Med.Statin.LLDyes 7.866e-01 1.271e+00 0.36892 1.6772
Med.all.antiplateletyes 9.045e-01 1.106e+00 0.29399 2.7825
GFR_MDRD 1.000e+00 9.997e-01 0.99977 1.0008
BMI 1.001e+00 9.989e-01 0.99990 1.0023
MedHx_CVDNo 5.534e-01 1.807e+00 0.25382 1.2066
stenose50-70% 6.126e+06 1.633e-07 0.00000 Inf
stenose70-90% 6.437e+06 1.553e-07 0.00000 Inf
stenose90-99% 6.875e+06 1.455e-07 0.00000 Inf
stenose100% (Occlusion) 8.516e-01 1.174e+00 0.00000 Inf
stenoseNA NA NA NA NA
stenose50-99% 1.737e-01 5.759e+00 0.00000 Inf
stenose70-99% NA NA NA NA
stenose99 NA NA NA NA
Concordance= 0.78 (se = 0.032 )
Likelihood ratio test= 37.58 on 31 df, p=0.2
Wald test = 26.88 on 31 df, p=0.7
Score (logrank) test = 36.13 on 31 df, p=0.2
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epstroke.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epstroke.3years
Protein...................: HDAC9
Effect size...............: -0.170258
Standard error............: 0.35101
Odds ratio (effect size)..: 0.843
Lower 95% CI..............: 0.424
Upper 95% CI..............: 1.678
T-value...................: -0.485054
P-value...................: 0.6276384
Sample size in model......: 530
Number of events..........: 37
* Analyzing the effect of plaque target-of-interest on [epcoronary.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
314 297
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 14,15,16,17,34,35 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 13,14,15,16,33,34 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = TEMP.DF)
n= 530, number of events= 42
(81 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -1.917e-01 8.256e-01 3.275e-01 -0.585 0.55831
Age 4.240e-03 1.004e+00 2.322e-02 0.183 0.85513
Gendermale 1.516e+00 4.554e+00 6.104e-01 2.483 0.01301 *
ORdate_year2003 -3.999e-01 6.704e-01 6.693e-01 -0.597 0.55018
ORdate_year2004 -1.227e+00 2.933e-01 7.828e-01 -1.567 0.11711
ORdate_year2005 -1.917e-01 8.256e-01 6.635e-01 -0.289 0.77266
ORdate_year2006 -3.560e-01 7.005e-01 6.501e-01 -0.548 0.58394
ORdate_year2007 -8.830e-01 4.135e-01 8.111e-01 -1.089 0.27628
ORdate_year2008 -4.915e-01 6.117e-01 7.423e-01 -0.662 0.50790
ORdate_year2009 -8.584e-01 4.238e-01 7.545e-01 -1.138 0.25525
ORdate_year2010 -1.852e+00 1.569e-01 1.148e+00 -1.613 0.10664
ORdate_year2011 -1.669e+00 1.885e-01 1.162e+00 -1.436 0.15109
ORdate_year2012 -9.168e-01 3.998e-01 1.167e+00 -0.786 0.43213
ORdate_year2013 -1.926e+01 4.303e-09 1.620e+04 -0.001 0.99905
ORdate_year2014 -1.794e+01 1.620e-08 3.501e+04 -0.001 0.99959
ORdate_year2015 -1.742e+01 2.726e-08 1.583e+04 -0.001 0.99912
ORdate_year2016 -1.742e+01 2.719e-08 3.501e+04 0.000 0.99960
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
ORdate_year2020 NA NA 0.000e+00 NA NA
ORdate_year2021 NA NA 0.000e+00 NA NA
ORdate_year2022 NA NA 0.000e+00 NA NA
Hypertension.compositeyes 1.573e+00 4.821e+00 1.031e+00 1.526 0.12702
DiabetesStatusDiabetes 3.327e-01 1.395e+00 3.492e-01 0.953 0.34068
SmokerStatusEx-smoker -2.464e-01 7.816e-01 3.558e-01 -0.693 0.48857
SmokerStatusNever smoked 1.276e-01 1.136e+00 4.770e-01 0.267 0.78908
Med.Statin.LLDyes 4.294e-01 1.536e+00 4.406e-01 0.974 0.32983
Med.all.antiplateletyes 2.228e-02 1.023e+00 5.714e-01 0.039 0.96890
GFR_MDRD -4.257e-04 9.996e-01 2.562e-04 -1.662 0.09661 .
BMI 4.625e-04 1.000e+00 5.521e-04 0.838 0.40220
MedHx_CVDNo -1.129e+00 3.234e-01 4.271e-01 -2.643 0.00822 **
stenose50-70% 1.097e+00 2.994e+00 1.877e+04 0.000 0.99995
stenose70-90% 1.896e+01 1.717e+08 1.797e+04 0.001 0.99916
stenose90-99% 1.884e+01 1.519e+08 1.797e+04 0.001 0.99916
stenose100% (Occlusion) 1.314e+00 3.722e+00 2.298e+04 0.000 0.99995
stenoseNA NA NA 0.000e+00 NA NA
stenose50-99% 1.578e+00 4.847e+00 3.935e+04 0.000 0.99997
stenose70-99% NA NA 0.000e+00 NA NA
stenose99 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 8.256e-01 1.211e+00 0.43451 1.569
Age 1.004e+00 9.958e-01 0.95957 1.051
Gendermale 4.554e+00 2.196e-01 1.37644 15.064
ORdate_year2003 6.704e-01 1.492e+00 0.18057 2.489
ORdate_year2004 2.933e-01 3.410e+00 0.06324 1.360
ORdate_year2005 8.256e-01 1.211e+00 0.22491 3.030
ORdate_year2006 7.005e-01 1.428e+00 0.19591 2.505
ORdate_year2007 4.135e-01 2.418e+00 0.08435 2.027
ORdate_year2008 6.117e-01 1.635e+00 0.14279 2.621
ORdate_year2009 4.238e-01 2.359e+00 0.09659 1.860
ORdate_year2010 1.569e-01 6.373e+00 0.01654 1.488
ORdate_year2011 1.885e-01 5.306e+00 0.01931 1.839
ORdate_year2012 3.998e-01 2.501e+00 0.04060 3.938
ORdate_year2013 4.303e-09 2.324e+08 0.00000 Inf
ORdate_year2014 1.620e-08 6.172e+07 0.00000 Inf
ORdate_year2015 2.726e-08 3.668e+07 0.00000 Inf
ORdate_year2016 2.719e-08 3.678e+07 0.00000 Inf
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
ORdate_year2020 NA NA NA NA
ORdate_year2021 NA NA NA NA
ORdate_year2022 NA NA NA NA
Hypertension.compositeyes 4.821e+00 2.074e-01 0.63929 36.356
DiabetesStatusDiabetes 1.395e+00 7.170e-01 0.70348 2.765
SmokerStatusEx-smoker 7.816e-01 1.279e+00 0.38919 1.570
SmokerStatusNever smoked 1.136e+00 8.802e-01 0.44605 2.894
Med.Statin.LLDyes 1.536e+00 6.509e-01 0.64775 3.644
Med.all.antiplateletyes 1.023e+00 9.780e-01 0.33366 3.134
GFR_MDRD 9.996e-01 1.000e+00 0.99907 1.000
BMI 1.000e+00 9.995e-01 0.99938 1.002
MedHx_CVDNo 3.234e-01 3.092e+00 0.14003 0.747
stenose50-70% 2.994e+00 3.340e-01 0.00000 Inf
stenose70-90% 1.717e+08 5.826e-09 0.00000 Inf
stenose90-99% 1.519e+08 6.583e-09 0.00000 Inf
stenose100% (Occlusion) 3.722e+00 2.687e-01 0.00000 Inf
stenoseNA NA NA NA NA
stenose50-99% 4.847e+00 2.063e-01 0.00000 Inf
stenose70-99% NA NA NA NA
stenose99 NA NA NA NA
Concordance= 0.784 (se = 0.026 )
Likelihood ratio test= 46.96 on 31 df, p=0.03
Wald test = 24.64 on 31 df, p=0.8
Score (logrank) test = 38.53 on 31 df, p=0.2
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epcoronary.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epcoronary.3years
Protein...................: HDAC9
Effect size...............: -0.191689
Standard error............: 0.327479
Odds ratio (effect size)..: 0.826
Lower 95% CI..............: 0.435
Upper 95% CI..............: 1.569
T-value...................: -0.585348
P-value...................: 0.5583135
Sample size in model......: 530
Number of events..........: 42
* Analyzing the effect of plaque target-of-interest on [epcvdeath.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
314 297
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 12,13,14,15,16,17,24,33,34,35,38 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 11,12,13,14,15,16,23,32,33,34,37 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = TEMP.DF)
n= 530, number of events= 23
(81 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -4.178e-01 6.585e-01 4.762e-01 -0.877 0.3803
Age 7.242e-02 1.075e+00 3.502e-02 2.068 0.0387 *
Gendermale 3.605e-01 1.434e+00 5.808e-01 0.621 0.5348
ORdate_year2003 -1.450e+00 2.346e-01 9.303e-01 -1.559 0.1191
ORdate_year2004 -1.161e+00 3.132e-01 9.034e-01 -1.285 0.1987
ORdate_year2005 -4.916e-01 6.117e-01 7.785e-01 -0.631 0.5278
ORdate_year2006 -1.264e+00 2.825e-01 9.076e-01 -1.393 0.1637
ORdate_year2007 -1.552e+00 2.118e-01 9.501e-01 -1.634 0.1023
ORdate_year2008 -5.706e-01 5.652e-01 8.332e-01 -0.685 0.4935
ORdate_year2009 -1.415e+00 2.430e-01 9.732e-01 -1.454 0.1460
ORdate_year2010 -5.791e-01 5.604e-01 9.477e-01 -0.611 0.5412
ORdate_year2011 -2.076e+01 9.612e-10 1.356e+04 -0.002 0.9988
ORdate_year2012 -2.036e+01 1.438e-09 2.044e+04 -0.001 0.9992
ORdate_year2013 -2.053e+01 1.209e-09 4.964e+04 0.000 0.9997
ORdate_year2014 -1.974e+01 2.686e-09 1.173e+05 0.000 0.9999
ORdate_year2015 -2.054e+01 1.202e-09 4.934e+04 0.000 0.9997
ORdate_year2016 -1.913e+01 4.896e-09 1.173e+05 0.000 0.9999
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
ORdate_year2020 NA NA 0.000e+00 NA NA
ORdate_year2021 NA NA 0.000e+00 NA NA
ORdate_year2022 NA NA 0.000e+00 NA NA
Hypertension.compositeyes 1.907e+01 1.918e+08 7.276e+03 0.003 0.9979
DiabetesStatusDiabetes 1.212e+00 3.361e+00 4.799e-01 2.526 0.0115 *
SmokerStatusEx-smoker -8.080e-01 4.457e-01 4.866e-01 -1.661 0.0968 .
SmokerStatusNever smoked -1.105e+00 3.313e-01 8.126e-01 -1.360 0.1740
Med.Statin.LLDyes -5.452e-01 5.797e-01 4.731e-01 -1.152 0.2492
Med.all.antiplateletyes -1.117e+00 3.272e-01 5.923e-01 -1.886 0.0593 .
GFR_MDRD -6.217e-04 9.994e-01 3.621e-04 -1.717 0.0860 .
BMI 8.061e-04 1.001e+00 8.100e-04 0.995 0.3197
MedHx_CVDNo -7.928e-01 4.526e-01 5.899e-01 -1.344 0.1790
stenose50-70% 1.913e+01 2.025e+08 3.916e+04 0.000 0.9996
stenose70-90% 1.905e+01 1.880e+08 3.916e+04 0.000 0.9996
stenose90-99% 1.914e+01 2.049e+08 3.916e+04 0.000 0.9996
stenose100% (Occlusion) 2.515e+00 1.237e+01 4.692e+04 0.000 1.0000
stenoseNA NA NA 0.000e+00 NA NA
stenose50-99% 3.679e+01 9.479e+15 1.256e+05 0.000 0.9998
stenose70-99% NA NA 0.000e+00 NA NA
stenose99 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 6.585e-01 1.519e+00 0.25897 1.674
Age 1.075e+00 9.301e-01 1.00379 1.151
Gendermale 1.434e+00 6.973e-01 0.45939 4.477
ORdate_year2003 2.346e-01 4.263e+00 0.03788 1.453
ORdate_year2004 3.132e-01 3.193e+00 0.05331 1.840
ORdate_year2005 6.117e-01 1.635e+00 0.13299 2.813
ORdate_year2006 2.825e-01 3.539e+00 0.04770 1.674
ORdate_year2007 2.118e-01 4.721e+00 0.03290 1.363
ORdate_year2008 5.652e-01 1.769e+00 0.11040 2.893
ORdate_year2009 2.430e-01 4.116e+00 0.03607 1.637
ORdate_year2010 5.604e-01 1.784e+00 0.08746 3.591
ORdate_year2011 9.612e-10 1.040e+09 0.00000 Inf
ORdate_year2012 1.438e-09 6.954e+08 0.00000 Inf
ORdate_year2013 1.209e-09 8.269e+08 0.00000 Inf
ORdate_year2014 2.686e-09 3.723e+08 0.00000 Inf
ORdate_year2015 1.202e-09 8.323e+08 0.00000 Inf
ORdate_year2016 4.896e-09 2.043e+08 0.00000 Inf
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
ORdate_year2020 NA NA NA NA
ORdate_year2021 NA NA NA NA
ORdate_year2022 NA NA NA NA
Hypertension.compositeyes 1.918e+08 5.214e-09 0.00000 Inf
DiabetesStatusDiabetes 3.361e+00 2.975e-01 1.31207 8.609
SmokerStatusEx-smoker 4.457e-01 2.243e+00 0.17174 1.157
SmokerStatusNever smoked 3.313e-01 3.018e+00 0.06738 1.629
Med.Statin.LLDyes 5.797e-01 1.725e+00 0.22937 1.465
Med.all.antiplateletyes 3.272e-01 3.056e+00 0.10249 1.045
GFR_MDRD 9.994e-01 1.001e+00 0.99867 1.000
BMI 1.001e+00 9.992e-01 0.99922 1.002
MedHx_CVDNo 4.526e-01 2.209e+00 0.14242 1.438
stenose50-70% 2.025e+08 4.938e-09 0.00000 Inf
stenose70-90% 1.880e+08 5.319e-09 0.00000 Inf
stenose90-99% 2.049e+08 4.880e-09 0.00000 Inf
stenose100% (Occlusion) 1.237e+01 8.084e-02 0.00000 Inf
stenoseNA NA NA NA NA
stenose50-99% 9.479e+15 1.055e-16 0.00000 Inf
stenose70-99% NA NA NA NA
stenose99 NA NA NA NA
Concordance= 0.862 (se = 0.031 )
Likelihood ratio test= 46.37 on 31 df, p=0.04
Wald test = 23.26 on 31 df, p=0.8
Score (logrank) test = 42.37 on 31 df, p=0.08
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epcvdeath.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epcvdeath.3years
Protein...................: HDAC9
Effect size...............: -0.417773
Standard error............: 0.476158
Odds ratio (effect size)..: 0.659
Lower 95% CI..............: 0.259
Upper 95% CI..............: 1.674
T-value...................: -0.877382
P-value...................: 0.3802791
Sample size in model......: 530
Number of events..........: 23 
cat("- Edit the column names...\n")- Edit the column names...colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")- Correct the variable types...COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
cat("- Writing results to Excel-file...\n")- Writing results to Excel-file...head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL2.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")- Removing intermediate files...rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times30)){
eptime = times30[i]
ep = endpoints30[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
ylim = c(0.75, 1),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [days]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".30days.pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 1 (Simple model)
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age+Gender + ORdate_year, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age+Gender + ORdate_year, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
ylim = c(0.75, 1), xlim = c(0,3), col = c("#595A5C", "#DB003F", "#1290D9"),
# ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [days]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL1.30days.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}
cat("- Edit the column names...\n")
colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")
COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
library(openxlsx)
cat("- Writing results to Excel-file...\n")
head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL1.30days.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colnames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")
rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times30)){
eptime = times30[i]
ep = endpoints30[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
#as.integer(TEMP.DF[,ep] == "Excluded")
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
ylim = c(0.75, 1),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [days]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".30days.pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 2 adjusted for age, sex, hypertension, diabetes, smoking, LDL-C levels, lipid-lowering drugs, antiplatelet drugs, eGFR, BMI, history of CVD, level of stenosis
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
ylim = c(0.75, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
# ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [days]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL2.30days.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}
cat("- Edit the column names...\n")
colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")
COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
cat("- Writing results to Excel-file...\n")
head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL2.30days.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")
rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times90)){
eptime = times90[i]
ep = endpoints90[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
ylim = c(0.75, 1),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [days]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".90days.pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 1 (Simple model)
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age+Gender + ORdate_year, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age+Gender + ORdate_year, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
ylim = c(0.75, 1), xlim = c(0,3), col = c("#595A5C", "#DB003F", "#1290D9"),
# ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [days]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL1.90days.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}
cat("- Edit the column names...\n")
colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")
COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
library(openxlsx)
cat("- Writing results to Excel-file...\n")
head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL1.90days.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")
rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times90)){
eptime = times90[i]
ep = endpoints90[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
#as.integer(TEMP.DF[,ep] == "Excluded")
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
ylim = c(0.75, 1),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [days]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".90days.pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 2 adjusted for age, sex, hypertension, diabetes, smoking, LDL-C levels, lipid-lowering drugs, antiplatelet drugs, eGFR, BMI, history of CVD, level of stenosis
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
ylim = c(0.75, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
# ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [days]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL2.90days.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}
cat("- Edit the column names...\n")
colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")
COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
cat("- Writing results to Excel-file...\n")
head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL2.90days.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")
rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)We correlated plaque levels of the biomarkers.
# Installation of ggcorrplot()
# --------------------------------
# if(!require(devtools))
# install.packages("devtools")
# devtools::install_github("kassambara/ggcorrplot")
library(ggcorrplot)Loading required package: ggplot2
# Creating matrix - inverse-rank transformation
# --------------------------------
# AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9,
# select = c("IL6_rank", "MCP1_rank", "IL6_pg_ug_2015_rank", "MCP1_pg_ug_2015_rank", "IL6R_pg_ug_2015_rank",
# TRAITS.BIN, TRAITS.CON.RANK)
# )
# AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9,
# select = c("MCP1_rank", "MCP1_pg_ug_2015_rank",
# TRAITS.BIN, TRAITS.CON.RANK)
# )
AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9,
select = c("HDAC9",
TRAITS.BIN,
TRAITS.CON.RANK,
"Symptoms.5G", "AsymptSympt", "EP_major", "EP_composite")
)
AERNASE.clin.hdac9.temp$CalcificationPlaque <- as.numeric(AERNASE.clin.hdac9.temp$CalcificationPlaque)
AERNASE.clin.hdac9.temp$CollagenPlaque <- as.numeric(AERNASE.clin.hdac9.temp$CollagenPlaque)
AERNASE.clin.hdac9.temp$Fat10Perc <- as.numeric(AERNASE.clin.hdac9.temp$Fat10Perc)
AERNASE.clin.hdac9.temp$MAC_binned <- as.numeric(AERNASE.clin.hdac9.temp$MAC_binned)
AERNASE.clin.hdac9.temp$SMC_binned <- as.numeric(AERNASE.clin.hdac9.temp$SMC_binned)
AERNASE.clin.hdac9.temp$IPH <- as.numeric(AERNASE.clin.hdac9.temp$IPH)
AERNASE.clin.hdac9.temp$Symptoms.5G <- as.numeric(AERNASE.clin.hdac9.temp$Symptoms.5G)
AERNASE.clin.hdac9.temp$AsymptSympt <- as.numeric(AERNASE.clin.hdac9.temp$AsymptSympt)
AERNASE.clin.hdac9.temp$EP_major <- as.numeric(AERNASE.clin.hdac9.temp$EP_major)
AERNASE.clin.hdac9.temp$EP_composite <- as.numeric(AERNASE.clin.hdac9.temp$EP_composite)
str(AERNASE.clin.hdac9.temp)'data.frame': 611 obs. of 15 variables:
$ HDAC9 : int 24 11 7 5 13 46 14 0 30 12 ...
$ CalcificationPlaque : num 2 2 2 1 1 1 2 2 2 1 ...
$ CollagenPlaque : num 2 2 2 2 2 2 1 2 2 2 ...
$ Fat10Perc : num 2 2 2 2 1 1 2 1 2 1 ...
$ IPH : num 1 2 2 1 1 1 1 2 1 1 ...
$ MAC_binned : num 2 1 2 1 1 2 2 1 2 1 ...
$ SMC_binned : num 2 2 2 2 2 2 1 2 2 2 ...
$ MAC_rankNorm : num -0.0697 -0.7511 0.9114 -0.38 -1.0968 ...
$ SMC_rankNorm : num 1.522 0.58 0.673 1.256 0.342 ...
$ MAC_SMC_ratio_rank : num -0.215 -0.88 0.728 -0.36 -1.31 ...
$ VesselDensity_rankNorm: num -0.723 0.769 0.854 1.14 -0.861 ...
$ Symptoms.5G : num 4 4 1 4 2 4 5 5 5 5 ...
$ AsymptSympt : num 2 2 1 2 1 2 2 2 2 2 ...
$ EP_major : num 1 1 1 1 2 2 2 1 1 1 ...
$ EP_composite : num 1 2 2 1 2 2 2 1 1 1 ...AERNASE.clin.hdac9.matrix.RANK <- as.matrix(AERNASE.clin.hdac9.temp)
rm(AERNASE.clin.hdac9.temp)corr_biomarkers.rank <- round(cor(AERNASE.clin.hdac9.matrix.RANK,
use = "pairwise.complete.obs", #the correlation or covariance between each pair of variables is computed using all complete pairs of observations on those variables
method = "spearman"), 3)
corr_biomarkers.rank HDAC9 CalcificationPlaque CollagenPlaque Fat10Perc IPH MAC_binned SMC_binned MAC_rankNorm SMC_rankNorm MAC_SMC_ratio_rank
HDAC9 1.000 0.044 -0.075 0.104 0.051 0.044 -0.095 -0.005 -0.153 -0.091
CalcificationPlaque 0.044 1.000 0.126 0.023 0.084 -0.079 0.057 -0.008 0.008 0.036
CollagenPlaque -0.075 0.126 1.000 -0.077 0.002 0.010 0.336 0.016 0.305 0.096
Fat10Perc 0.104 0.023 -0.077 1.000 0.204 0.294 -0.157 0.196 -0.253 -0.014
IPH 0.051 0.084 0.002 0.204 1.000 0.094 -0.118 0.135 -0.086 -0.025
MAC_binned 0.044 -0.079 0.010 0.294 0.094 1.000 0.043 0.628 -0.040 -0.052
SMC_binned -0.095 0.057 0.336 -0.157 -0.118 0.043 1.000 0.110 0.557 0.083
MAC_rankNorm -0.005 -0.008 0.016 0.196 0.135 0.628 0.110 1.000 0.213 -0.027
SMC_rankNorm -0.153 0.008 0.305 -0.253 -0.086 -0.040 0.557 0.213 1.000 0.176
MAC_SMC_ratio_rank -0.091 0.036 0.096 -0.014 -0.025 -0.052 0.083 -0.027 0.176 1.000
VesselDensity_rankNorm -0.031 -0.045 0.169 0.023 0.200 0.077 0.115 0.088 0.083 0.023
Symptoms.5G -0.070 0.027 -0.029 0.059 0.027 -0.011 -0.018 0.015 0.015 -0.004
AsymptSympt -0.024 -0.025 -0.077 0.126 0.044 0.017 -0.097 0.029 -0.067 -0.005
EP_major -0.062 0.059 0.045 -0.021 0.006 -0.029 0.077 -0.002 0.032 -0.009
EP_composite -0.029 0.066 0.040 -0.041 0.071 -0.063 0.028 -0.014 0.028 0.032
VesselDensity_rankNorm Symptoms.5G AsymptSympt EP_major EP_composite
HDAC9 -0.031 -0.070 -0.024 -0.062 -0.029
CalcificationPlaque -0.045 0.027 -0.025 0.059 0.066
CollagenPlaque 0.169 -0.029 -0.077 0.045 0.040
Fat10Perc 0.023 0.059 0.126 -0.021 -0.041
IPH 0.200 0.027 0.044 0.006 0.071
MAC_binned 0.077 -0.011 0.017 -0.029 -0.063
SMC_binned 0.115 -0.018 -0.097 0.077 0.028
MAC_rankNorm 0.088 0.015 0.029 -0.002 -0.014
SMC_rankNorm 0.083 0.015 -0.067 0.032 0.028
MAC_SMC_ratio_rank 0.023 -0.004 -0.005 -0.009 0.032
VesselDensity_rankNorm 1.000 -0.088 -0.101 0.119 0.171
Symptoms.5G -0.088 1.000 0.884 0.018 0.053
AsymptSympt -0.101 0.884 1.000 0.016 0.036
EP_major 0.119 0.018 0.016 1.000 0.685
EP_composite 0.171 0.053 0.036 0.685 1.000corr_biomarkers_p.rank <- ggcorrplot::cor_pmat(AERNASE.clin.hdac9.matrix.RANK, use = "pairwise.complete.obs", method = "spearman", exact = FALSE)# Add correlation coefficients
# --------------------------------
# argument lab = TRUE
ggcorrplot(corr_biomarkers.rank,
method = "square",
type = "lower",
title = "Cross biomarker correlations",
show.legend = TRUE, legend.title = bquote("Spearman's"~italic(rho)),
ggtheme = ggplot2::theme_minimal, outline.color = "#FFFFFF",
show.diag = TRUE,
hc.order = FALSE,
lab = FALSE,
digits = 3,
# p.mat = corr_biomarkers_p.rank, sig.level = 0.05,
colors = c("#1290D9", "#FFFFFF", "#E55738"))
library(data.table)Registered S3 method overwritten by 'data.table':
method from
print.data.table
data.table 1.14.2 using 1 threads (see ?getDTthreads). Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********
Attaching package: ‘data.table’
The following object is masked from ‘package:SummarizedExperiment’:
shift
The following object is masked from ‘package:GenomicRanges’:
shift
The following object is masked from ‘package:IRanges’:
shift
The following objects are masked from ‘package:S4Vectors’:
first, second# flattenCorrMatrix
# --------------------------------
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
ut <- upper.tri(cormat)
data.frame(
biomarker_row = rownames(cormat)[row(cormat)[ut]],
biomarker_column = rownames(cormat)[col(cormat)[ut]],
spearman_cor =(cormat)[ut],
pval = pmat[ut]
)
}
corr_biomarkers.rank.df <- as.data.table(flattenCorrMatrix(corr_biomarkers.rank, corr_biomarkers_p.rank))
DT::datatable(corr_biomarkers.rank.df)NA# chart of a correlation matrix
# --------------------------------
# Alternative solution https://www.r-graph-gallery.com/199-correlation-matrix-with-ggally.html
install.packages.auto("PerformanceAnalytics")Installing package into ‘/Users/swvanderlaan/Library/R/x86_64/4.1/library’
(as ‘lib’ is unspecified)
also installing the dependencies ‘RPEIF’, ‘RPEGLMEN’, ‘pyinit’, ‘googleVis’, ‘RPESE’, ‘RobStatTM’
trying URL 'https://cloud.r-project.org/bin/macosx/contrib/4.1/RPEIF_1.2.4.tgz'
Content type 'application/x-gzip' length 487987 bytes (476 KB)
==================================================
downloaded 476 KB
trying URL 'https://cloud.r-project.org/bin/macosx/contrib/4.1/RPEGLMEN_1.1.1.tgz'
Content type 'application/x-gzip' length 1278219 bytes (1.2 MB)
==================================================
downloaded 1.2 MB
trying URL 'https://cloud.r-project.org/bin/macosx/contrib/4.1/pyinit_1.1.2.tgz'
Content type 'application/x-gzip' length 53583 bytes (52 KB)
==================================================
downloaded 52 KB
trying URL 'https://cloud.r-project.org/bin/macosx/contrib/4.1/googleVis_0.6.11.tgz'
Content type 'application/x-gzip' length 483799 bytes (472 KB)
==================================================
downloaded 472 KB
trying URL 'https://cloud.r-project.org/bin/macosx/contrib/4.1/RPESE_1.2.4.tgz'
Content type 'application/x-gzip' length 340165 bytes (332 KB)
==================================================
downloaded 332 KB
trying URL 'https://cloud.r-project.org/bin/macosx/contrib/4.1/RobStatTM_1.0.3.tgz'
Content type 'application/x-gzip' length 778852 bytes (760 KB)
==================================================
downloaded 760 KB
trying URL 'https://cloud.r-project.org/bin/macosx/contrib/4.1/PerformanceAnalytics_2.0.4.tgz'
Content type 'application/x-gzip' length 3155393 bytes (3.0 MB)
==================================================
downloaded 3.0 MB
The downloaded binary packages are in
/var/folders/cj/1vxt4xb11m1cx7wn020f8hww0000gn/T//Rtmpj0vR9W/downloaded_packagesLoading required package: PerformanceAnalytics
Loading required package: xts
Loading required package: zoo
Attaching package: ‘zoo’
The following objects are masked from ‘package:base’:
as.Date, as.Date.numeric
Attaching package: ‘xts’
The following objects are masked from ‘package:data.table’:
first, last
The following object is masked from ‘package:S4Vectors’:
first
Attaching package: ‘PerformanceAnalytics’
The following object is masked from ‘package:graphics’:
legendchart.Correlation.new <- function (R, histogram = TRUE, method = c("pearson", "kendall",
"spearman"), ...)
{
x = checkData(R, method = "matrix")
if (missing(method))
method = method[1]
cormeth <- method
panel.cor <- function(x, y, digits = 2, prefix = "", use = "pairwise.complete.obs",
method = cormeth, cex.cor, ...) {
usr <- par("usr")
on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- cor(x, y, use = use, method = method)
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste(prefix, txt, sep = "")
if (missing(cex.cor))
cex <- 0.8/strwidth(txt)
test <- cor.test(as.numeric(x), as.numeric(y), method = method)
Signif <- symnum(test$p.value, corr = FALSE, na = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), symbols = c("***",
"**", "*", ".", " "))
text(0.5, 0.5, txt, cex = cex * (abs(r) + 0.3)/1.3)
text(0.8, 0.8, Signif, cex = cex, col = 2)
}
f <- function(t) {
dnorm(t, mean = mean(x), sd = sd.xts(x))
}
dotargs <- list(...)
dotargs$method <- NULL
rm(method)
hist.panel = function(x, ... = NULL) {
par(new = TRUE)
hist(x, col = "#1290D9", probability = TRUE, axes = FALSE,
# hist(x, col = "light gray", probability = TRUE, axes = FALSE,
main = "", breaks = "FD")
lines(density(x, na.rm = TRUE), col = "#E55738", lwd = 1)
rug(x)
}
if (histogram)
pairs(x, gap = 0, lower.panel = panel.smooth, upper.panel = panel.cor,
diag.panel = hist.panel, ...)
else pairs(x, gap = 0, lower.panel = panel.smooth, upper.panel = panel.cor, ...)
}
chart.Correlation.new(AERNASE.clin.hdac9.matrix.RANK, method = "spearman", histogram = TRUE, pch = 3)
# alternative chart of a correlation matrix
# --------------------------------
# Alternative solution https://www.r-graph-gallery.com/199-correlation-matrix-with-ggally.html
install.packages.auto("GGally")
# Quick display of two cabapilities of GGally, to assess the distribution and correlation of variables
library(GGally)
# From the help page:
ggpairs(AERNASE.clin.hdac9,
columns = c("HDAC9", TRAITS.BIN, TRAITS.CON.RANK, "Symptoms.5G", "AsymptSympt", "EP_major", "EP_composite"),
columnLabels = c("HDAC9",
"Calcification", "Collagen", "Fat 10%", "IPH", "Macrophages (binned)", "SMC (binned)", "Macrophages", "SMC", "Macrophage/SMC", "Vessel density",
"Symptoms", "Symptoms (grouped)", "MACE", "Composite"),
method = c("spearman"),
# ggplot2::aes(colour = Gender),
progress = FALSE)`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Version: v1.0.0
Last update: 2021-03-18
Written by: Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
Description: Script to analyse HDAC9 from the Ather-Express Biobank Study.
Minimum requirements: R version 3.5.2 (2018-12-20) -- 'Eggshell Igloo', macOS Mojave (10.14.2).
**MoSCoW To-Do List**
The things we Must, Should, Could, and Would have given the time we have.
_M_
_S_
_C_
_W_
**Changes log**
* v1.0.0 Inital version.sessionInfo()R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Monterey 12.2.1
Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] GGally_2.1.2 PerformanceAnalytics_2.0.4 xts_0.12.1 zoo_1.8-9 data.table_1.14.2
[6] ggcorrplot_0.1.3.999 ggplot2_3.3.5 rmarkdown_2.11 worcs_0.1.9.1 SummarizedExperiment_1.24.0
[11] Biobase_2.54.0 GenomicRanges_1.46.1 GenomeInfoDb_1.30.0 IRanges_2.28.0 S4Vectors_0.32.3
[16] BiocGenerics_0.40.0 MatrixGenerics_1.6.0 matrixStats_0.61.0
loaded via a namespace (and not attached):
[1] backports_1.4.1 BiocFileCache_2.2.0 plyr_1.8.6 lazyeval_0.2.2 BiocParallel_1.28.3 crosstalk_1.2.0
[7] usethis_2.1.5 digest_0.6.29 ensembldb_2.18.2 htmltools_0.5.2 fansi_1.0.2 magrittr_2.0.2
[13] memoise_2.0.1 gert_1.5.0 remotes_2.4.2 credentials_1.3.2 Biostrings_2.62.0 askpass_1.1
[19] prettyunits_1.1.1 colorspace_2.0-3 rappdirs_0.3.3 blob_1.2.2 xfun_0.29 dplyr_1.0.7
[25] callr_3.7.0 crayon_1.5.0 RCurl_1.98-1.5 jsonlite_1.7.2 glue_1.6.2 prereg_0.5.0
[31] gtable_0.3.0 zlibbioc_1.40.0 XVector_0.34.0 DelayedArray_0.20.0 car_3.0-12 pkgbuild_1.3.1
[37] abind_1.4-5 scales_1.1.1 DBI_1.1.2 rstatix_0.7.0 Rcpp_1.0.8.3 progress_1.2.2
[43] bit_4.0.4 DT_0.20 htmlwidgets_1.5.4 httr_1.4.2 RColorBrewer_1.1-2 ellipsis_0.3.2
[49] pkgconfig_2.0.3 reshape_0.8.8 XML_3.99-0.8 farver_2.1.0 sass_0.4.0 dbplyr_2.1.1
[55] utf8_1.2.2 tidyselect_1.1.1 labeling_0.4.2 rlang_1.0.2 reshape2_1.4.4 AnnotationDbi_1.56.2
[61] munsell_0.5.0 tools_4.1.2 cachem_1.0.6 cli_3.2.0 generics_0.1.1 RSQLite_2.2.9
[67] ranger_0.13.1 devtools_2.4.3 broom_0.7.11 evaluate_0.14 stringr_1.4.0 fastmap_1.1.0
[73] yaml_2.2.1 sys_3.4 rticles_0.22 processx_3.5.2 knitr_1.37 bit64_4.0.5
[79] fs_1.5.2 purrr_0.3.4 AnnotationFilter_1.18.0 KEGGREST_1.34.0 xml2_1.3.3 biomaRt_2.50.2
[85] compiler_4.1.2 rstudioapi_0.13 filelock_1.0.2 curl_4.3.2 png_0.1-7 testthat_3.1.1
[91] ggsignif_0.6.3 tibble_3.1.6 bslib_0.3.1 stringi_1.7.6 ps_1.6.0 method
[97] desc_1.4.0 lattice_0.20-45 ProtGenerics_1.26.0 Matrix_1.4-0 vctrs_0.3.8 pillar_1.7.0
[103] lifecycle_1.0.1 jquerylib_0.1.4 bitops_1.0-7 rtracklayer_1.54.0 R6_2.5.1 BiocIO_1.4.0
[109] sessioninfo_1.2.2 assertthat_0.2.1 pkgload_1.2.4 openssl_1.4.6 rprojroot_2.0.2 rjson_0.2.21
[115] withr_2.5.0 GenomicAlignments_1.30.0 Rsamtools_2.10.0 GenomeInfoDbData_1.2.7 parallel_4.1.2 hms_1.1.1
[121] quadprog_1.5-8 grid_4.1.2 tidyr_1.1.4 carData_3.0-5 ggpubr_0.4.0 restfulr_0.0.13 save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".bulkRNAseq.main_analysis.RData"))| © 1979-2022 Sander W. van der Laan | s.w.vanderlaan[at]gmail.com swvanderlaan.github.io. |